Pan Liu

Sidy Fall

Olivier Baledent

Accuracy of Real-Time Echo-Planar Imaging Phase Contrast MRI

Compared with CINE phase contrast MRI (CINE-PC), echo-planar imaging phase contrast (EPI-PC) can achieve realtime quantification of blood flow, with lower SNR. In this study, the pulsating real model of the simulated cerebral vasculature was used to verify the accuracy of EPI-PC. The imaging time of EPI-PC was 62ms/image at 100*60 spatial resolution. The reconstructed EPI-PC flow curve was extracted by homemade post-processing software. After comparison with the CINE-PC flow curve, it was concluded that EPI-PC can provide an average flow with less than 3% error, and its flow curve will be similar to the CINE-PC flow curve in shape.

Introduction

Today, CINE phase-contrast MRI (CINE-PC) [START_REF] Nj | Phase Contrast Cine Magnetic Resonance Imaging[END_REF] is the primary noninvasive method for quantifying blood flow and cerebrospinal fluid (CSF) in the cerebrovascular system by using cardiac synchronization [START_REF] Balédent | Cerebrospinal fluid dynamics and relation with blood flow: a magnetic resonance study with semiautomated cerebrospinal fluid segmentation[END_REF]. However, it cannot provide the effect of respiration on blood/LCS flow [START_REF] Chen | Dynamics of Respiratory and Cardiac CSF Motion Revealed with Real-Time Simultaneous Multi-Slice EPI Velocity Phase Contrast Imaging[END_REF] as it can only provide the average cardiac cycle reconstructed using all cardiac cycles obtained during the acquisition. Echo-planar imaging phase contrast (EPI-PC) [START_REF] Pan | Post Processing Software for Echo Planar Imaging Phase Contrast Sequence[END_REF] is a recent technique can achieve real-time quantification of blood/LCS flow and does not require synchronization. Meanwhile, EPI-PC is sensitive to noise and has a low signal-to-noise ratio (SNR). The purpose of this study was to validate the accuracy of EPI-PC in vitro and to quantify the effect of pixel size on both sequences.

Material

The study was performed on a clinical 3T scanner (Philips Achieva, maximum gradient 80mT/m, gradient increase rate 120mT/m/ms) with a 32-channel head coil.

A flow phantom consisting of 4 tubes in series with different diameters (Φ:9.5/6.4/4.4/2 mm) was attached to the pump to mimicking the vessels in craniospinal system (Figure .1), and a static water-filled tube was placed next to the tubes for calibration. The water pump provides pulsatile flow at 99 bpm with an average flow rate of 1150 mm 3 /s. This experiment will image the first tube (Φ9.5 mm) with CINE-PC and EPI-PC. CINE-PC provides 32 phase images to show an average pulse cycle. And EPI-PC contains 150 phase images to show the flow rate in real time (Figure 2).

Post-processing is performed by home-made software [START_REF] Chao | Accuracy of Phase-Contrast Flow Measurements in the Presence of Partial-Volume Effects[END_REF] to extract the flow curves. Calibration is performed by measuring the average velocity of the static tube. In addition, for comparison with the CINE-PC, The EPI-PC flow signal is reconstructed as flow curve with 32 sampling points in the same format as the CINE-PC flow curve.

METHODS

Both sequences were imaged 10 times using the default parameters (Figure 3) for tube-1. The average flow rate of the pump (1150 mm 3 /s) and the cross-sectional area of the tube (70.8 mm 2 ) were used as gold standards, and the confidence intervals for the average flow rate and the segmentation area were set to 100 ± 10% (average flow: [1035 mm 3 /s: 1265 mm 3 /s]; segmentation area: [63.72 mm 2 : 77.8 mm 2 ]). The Bland & Altman test was then used to verify the similarity of the shapes of these two curves.

Based on the default parameters, pixel size was used as a variable for 10 group images (from 0.8 mm to 4.4 mm with 0.4 mm interval), and each group was repeated four times to observe its effect on the two sequences.

Results

After 10 measurements, reconstructed EPI-PC flow curves were synchronized with CINE-PC flow curves at the peak flow, as shown in Figure 4 A. The average flow of EPI-PC was 1116 ± 24.5 mm 3 /s with a coefficient of variation of 2.3%. The average flow of CINE-PC was 1239 ± 26.3 mm3/s with a coefficient of variation of 2.3%. Both the average flow and the segmentation area (EPI-PC is 70.1 mm 2 and CINE-PC is 70 mm 2 ) are within the confidence interval. The Bland-Altman plot (Figure 4 B) reflects that the differences of the 32 points corresponding to the two curves are within the limits of agreement (LoA), which can indicate that the two curves have similar shapes.

Figure 5 shows that the average flow of EPI-PC is within the confidence interval when the pixel size is between 1.2 mm and 2.8 mm, while the pixel size of CINE-PC needs to be set between 1.2 mm and 2.4 mm to provide an accurate average flow.

Discussion

Under the present experimental conditions, the reconstructed EPI-PC flow curve can replace the CINE-PC flow curve. However, in vivo, cerebral arterial blood flow waves usually have higher frequency components (4Hz-5Hz), and for the current sampling frequency, EPI-PC may have difficulty in correctly representing these high frequency waveforms. CINE-PC completes multiple phase encodings in each pulse cycle and fills them in K-space of the corresponding phase images, so that even if the acquisition time increases, it does not affect its pseudo-sampling interval ∆t. However, for EPI-PC with real-time imaging, ∆t is directly related to the acquisition time. Therefore, increasing the sampling frequency can improve the accuracy of the reconstructed curve shape without affecting the accuracy of the results. In the case of larger pixel size, not only the signal-to-noise ratio can be improved, but also the sampling frequency can be increased. However, if the pixel size is too large, a partial volume effect will occur, which distorts the results. The findings of this study indicate that setting the pixel size of EPI-PC between 25%-30% of the target vessel/LCS diameter can obtain a higher signal-to-noise ratio and a higher temporal resolution.

Conclusion

This experiment demonstrates the feasibility of using EPI-PC for in vivo real-time imaging. After post-processing with specific software, EPI-PC can provide the same flow curve as CINE-PC and obtain time-dimensional data. Increasing the sampling frequency allows for better display of the details of the reconstructed flow curve from EPI-PC. With the application of more acceleration techniques such as compressed sensing, EPI sequences will have greater potential for clinical diagnosis and for studying the effect of respiration on blood flow.

Figures

  Figures

Figure 1 :

 1 Figure 1: Left: A pulsatile flow phantom simulating the vasculature of the cranial system is located in the center of the head coil. An elastic balloon is connected in parallel to the outflow tube for connection to the synchronizer of the CINE-PC. Right: Amplitude and phase images of the 4 tubes and the static water-filled tube in the imaging plane.

Figure 2 :

 2 Figure2: CINE-PC phase image A) and its calibrated flow curve B). EPI-PC phase image C) and its calibrated flow curve D) with the minimum points (red cycles) in each cycle automatically found by the software which will be used as the 'cutting' position, all the pulse cycles of EPI-PC will be used to reconstruct an average pulsed cycle E).

Figure 3 :Figure 4 :

 34 Figure 3: Default parameters of CINE-PC and EPI-PC

Figure 5 :

 5 Figure 5: Results (segmentation area & average flow) distribution of EPI-PC (red points) and CINE-PC (blue points) under different pixel size (different color saturation). The X axis represents the segmentation area, the purple shade is the confidence interval of segmentation area. The Y axis represent the flow rate, the green shade is the confidence interval of flow rate.

Acknowledgements

This research was supported by EquipEX FIGURES (Facing Faces Institute Guilding Research), European Union Interreg REVERT Project, Hanuman ANR-18-CE45-0014 and Region Haut de France.

Thanks to the staff members at the Facing Faces Institute (Amiens, France) for technical assistance.

Thanks to David Chechin from Phillips industry for his scientific support.


Félix Bigand
email: felix.bigand@lisn.upsaclay.fr

Elise Prigent
email: elise.prigent@lisn.upsaclay.fr

Annelies Braffort
email: annelies.braffort@lisn.upsaclay.fr

Synthesis for the Kinematic Control of Identity in Sign Language

Keywords: Sign Language, Anonymized Content, Identity Conversion, Motion Generation, Machine Learning

Sign Language (SL) animations generated from motion capture (mocap) of real signers convey critical information about their identity. It has been suggested that this information is mostly carried by statistics of the movements kinematics. Manipulating these statistics in the generation of SL movements could allow controlling the identity of the signer, notably to preserve anonymity. This paper tests this hypothesis by presenting a novel synthesis algorithm that manipulates the identity-specific statistics of mocap recordings. The algorithm produced convincing new versions of French Sign Language discourses, which accurately modulated the identity prediction of a machine learning model. These results open up promising perspectives toward the automatic control of identity in the motion animation of virtual signers.

Introduction

Using motion capture (mocap) systems, the movements of signers can be recorded with high accuracy and be used to produce natural and comprehensible content [START_REF] Lu | Collecting a motion-capture corpus of american sign language for data-driven generation research[END_REF][START_REF] Gibet | Building french sign language motion capture corpora for signing avatars[END_REF]. However, this process raises an unexpected problem, related to the human ability to identify individuals from their movements [START_REF] Troje | Person identification from biological motion: Effects of structural and kinematic cues[END_REF][START_REF] Loula | Recognizing people from their movement[END_REF][START_REF] Bläsing | My action, my self: Recognition of self-created but visually unfamiliar dance-like actions from point-light displays[END_REF]. As for spoken languages in the auditory domain, where voice parameters inform about the identity of a speaker, signers can be identified from their movements [START_REF] Bigand | Person identification based on sign language motion: Insights from human perception and computational modeling[END_REF]. We present a synthesis algorithm for controlling the motion features that characterize the identity of a signer. This would allow producing anonymized, non-identifiable, content with virtual signers, which is crucial (e.g., for sharing anonymized testimony) given that Sign Languages (SLs) have no written form [START_REF] Lee | American sign language video anonymization to support online participation of deaf and hard of hearing users[END_REF].

In line with prior work on non-SL movements (Troje et al., 2005;[START_REF] Carlson | Dance to your own drum: Identification of musical genre and individual dancer from motion capture using machine learning[END_REF][START_REF] Zhang | View-independent person identification from human gait[END_REF], our recent studies suggested that identity was mainly inferred from the kinematic aspects of the movements, beyond size, shape or posture of the signers [START_REF] Bigand | Person identification based on sign language motion: Insights from human perception and computational modeling[END_REF][START_REF] Bigand | Machine learning of motion statistics reveals the kinematic signature of a person's identity in sign language[END_REF]. Using a machine learning model, we automatically extracted the specific kinematic aspects of motion that carry identity, using time-averaged statistics (Section 2). The present synthesis algorithm was then developed in order to manipulate the identity-specific statistics of original mocap recordings (Section 3). We tested the performance of the synthesis algorithm by modifying the identity attribute of mocap recordings in French Sign Language, and by assessing the identity inferred from the new excerpts (Section 4). This constitutes the first step toward automatically anonymizing the movements of signers in SL animations, in the same way as for the voice of a speaker, which can be anonymized by modifying specific vocal parameters (Section 5).

Motion statistics of identity

Mocap recordings were taken from MOCAP1 corpus [START_REF] Benchiheub | MOCAP1 corpus[END_REF]. Six signers had freely described the content of 24 pictures using French Sign Language (LSF). From each of the 24 original recordings, one mocap recording unit of 5-seconds duration was extracted from the beginning of the utterance (see examples in Videos 7.1 to 7.6). As shown in Figure 1 Statistics of the mocap examples were then computed as follows. Based on previous research investigating the perception of auditory and visual textures [START_REF] Mcdermott | Sound texture perception via statistics of the auditory periphery: evidence from sound synthesis[END_REF][START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF], we measured the first two moments (i.e., mean and standard deviation (SD)) of position and velocity, and covariances of velocity between body markers. The first two moments of position and velocity described their statistical distributions, which may vary from one individual to another, as shown for expert gesture analysis [START_REF] Tits | Expert Gesture Analysis through Motion Capture using Statistical Modeling and Machine Learning[END_REF]. Moreover, the covariance of velocity allowed for quantifying the extent to which any two markers covaried with each other, in two directions. This latter statistic has been shown to allow for automatic person identification from dance movements [START_REF] Carlson | Dance to your own drum: Identification of musical genre and individual dancer from motion capture using machine learning[END_REF]. These statistics vary substantially across the mocap data of different signers. For instance, as shown in Figure 2, the position and velocity data of one body marker are distributed differently across signers, for one mocap example (i.e., for comparable content: the description of the same picture in LSF). Distributions of position data differ in location of the peak (captured by the mean) and width (captured by the standard deviation). Figure 3 further supports that the two moments of position and velocity may capture substantial differences across signers. Furthermore, velocity covariances capture different aspects of motor coordination between the markers in three dimensions, which can differ across signers. Various distinct coordination patterns can be extracted. For instance, for mocap example 24, the movements of Signer 2 show an overall substantial (positive) covariance between body markers along the Y axis, while this covariance is near zero for Signer 4 (Figure 4.A). Inversely, Signer 4 displays an important (negative) covariance of movements of the right arm and hand along the Y axis with the trunk (i.e., stomach and sternum) and head markers along the X axis, while this covariance is less important for Signer 2 (Figure 4.B). Taken together, these examples raise the possibility that the identity of a signer is conveyed by statistical properties of his or her movements.

Methods

The automatic signer identification model presented in [START_REF] Bigand | Machine learning of motion statistics reveals the kinematic signature of a person's identity in sign language[END_REF] allowed extracting specific kinematic statistics that carry identity information about the signers. A linear classifier was trained to extract the statistics of the mocap data characteristic of identity (i.e., the ones that allow for accurate signer identification). Then, the aim of the present synthesis algorithm was to manipulate the statistics of an original SL mocap recording (i.e., impose new statistics to the original recording), in order to reduce (α < 0) or exaggerate (α > 0) the identity attribute, following Equation 1:

dα = d orig + αd k (1)
where dα is a vector containing the new target statistics to be imposed by the synthesis alogrithm, d orig is a vector containing the original statistics of the mocap example, d k is a vector containing the overall statistical patterns characteristic of the identity of Signer k, and α is a scalar related to the amount of reduction (α < 0) or exaggeration (α > 0) of the identity attribute.

The different steps of the synthesis process are displayed in Figure 5. In summary, the synthesis process consisted of modifying (i.e., "re-synthezing") an existing mocap recording in order to change the identity attribute of the signer, according to the following steps. First, statistics of the original mocap example are measured, while the discriminant statistical kinematic patterns are extracted by the automatic identification model (see [START_REF] Bigand | Machine learning of motion statistics reveals the kinematic signature of a person's identity in sign language[END_REF]). Then, the discriminant statistics characteristic of Signer k are either added to (α > 0) or subtracted from (α < 0) the ones of the original example (see Equation 1). Multiple manipulations can then be done using this technique, depending on the values of k and α. For instance, if the original mocap example relates to Signer 1, reducing the importance of her identity-specific statistics (i.e., k = 1, α < 0) would make her less identifiable (i.e., kinematic anonymization). By contrast, increasing the importance of the identity-specific statistics of Signer 2 (i.e., k = 2, α > 0) would make this latter signer identifiable while the SL movements were originally executed by Signer 1 (i.e., kinematic identity conversion).

Once the target statistics defined, they are imposed to the original mocap signal by the algorithm, which creates a new mocap excerpt.

Target statistics were imposed using an iterative process where a synthesized mocap signal (initialized with the content of the original mocap recording) is modified until its statistics are sufficiently close to the target ones dα . Mathematically, the objective of this process is to minimize the loss function that calculates the mean square of the differences between the target statistics and the statistics of the synthesized movements (see Equation 2). We imposed the first two moments (mean and SD) of position and velocity data and the covariance of velocity between markers, as they were found to be the most important statistics for signer identification [START_REF] Bigand | Extracting human characteristics from motion using machine learning : the case of identity in Sign Language[END_REF]. Imposing the mean of position and mean of velocity of the markers was done to maintain consistent motion data when synthesizing (e.g., to avoid the generation of unrealistic, non-biological, movements), although these two statistics had only minor role in the identification.

(C i,j,targ -C i,j,synth ) 2 loss tot = 5 i=1 loss i (2)
where µ pos,m , σ pos,m , µ vel,m and σ vel,m are the first two moments of position and velocity data of marker m (m ∈ [1, 54]), C i,j is the covariance of velocity between markers i and j. targ and synth subscripts distinguish between target statistics and statistics of the synthesized movements, respectively. In order to be able to minimize all of the five loss components of Equation 2 despite the differences in ranges of amplitude across statistics, we used a weighted loss function, whose weights then need to be optimized (see Equation 3). The loss function was then minimized using the Adam optimization algorithm for gradient descent. Each iterative step of the gradient descent modified the synthesized mocap signals (i.e., position temporal curves of the 19 markers along the three dimensions) so that they approached the target statistics.

loss tot = 5 i=1 w i loss i (3)
where loss i is the loss function related to one statistical measure and w i the optimized weight.

Initially, there was no constraint in the synthesis process that forced the position and velocity signals of the synthesized movements to remain consistent with their initial temporal structure in the original movements. The limitation of this first version of the algorithm is that, although it managed to impose the statistics present in Equation 2, the modifications applied to the new movements seemed to generate noise artifacts rather than changing relevant aspects of the motion of the signer (see Video 10.1). In fact, the imposing algorithm managed to impose the target statistics but by modifying the movements in an undesired manner. First, low-energy segments of the motion were modified in the same way as high-energy ones, which is not relevant as they may not be perceived by observers. Moreover, reaching the target statistics caused very rapid oscillations in the synthesized velocity temporal curves, which are unlikely to be perceived as biological motion by the observers (but rather noisy, wobbling, markers).

In order to modify the movements in proportion to their energy (i.e., modify the aspects of the movement at relevant times of actual, perceptible, motion), we included another target statistic in the imposing algorithm: the correlation of velocity between the original and synthesized movements. The algorithm then aimed to minimize the mean squared error between this correlation and a value of 1, which characterizes two signals that are perfectly positively correlated (see Equation 4). In other words, imposing this additional statistic (Equation 5) allowed forcing the velocity curves of the synthesized movements to be consistent with their initial temporal structure in the original mocap recording (Figure 6).

loss 6 = m (1 -ρ vel,m,synth ) 2 (4)
loss tot = 6 i=1 w i loss i (5)
where ρ vel,m is the correlation of velocity between the original and synthesized movements of marker m (m ∈ [1, 54]). The target correlation value is set to 1 for all markers, in order to preserve the original temporal structure of velocity curves.

Results

This synthesis procedure was run on mocap examples of different signers and for different modifications of the identity attribute. In order to visualize how these For instance, the movements of Signer 1 were modified so that the perceived identity was that of Signer 2 (i.e., identity conversion) (see Videos 10.3 and 10.4). Then, they were modified to make Signer 1 not identifiable, without making another signer identifiable specifically (i.e., anonymization) (see Videos 10.5 and 10.6). One further synthesis example of identity conversion (from Signer 2 to Signer 1) can be found in [START_REF] Bigand | Extracting human characteristics from motion using machine learning : the case of identity in Sign Language[END_REF].

In order to assess the extent to which the novel movements generated by our algorithm could convey a modified identity attribute (e.g., could be anonymized, or identified as movements of another signer), we tested our automatic signer identification model on the synthesized mocap examples. If the identity-specific aspects of the movements are correctly modified by the synthesis algorithm, then automatic identification from these synthesized examples should be compromised.

When converting the identity of Signer 1 into that of Signer 2, the automatic signer identification model identified the synthesized mocap example as that of Signer 2, while it identified the original motion as produced by Signer 1 (see Table 1). Then, when anonymizing the content of Signer 1, the signer identification model did not manage to identify Signer 1 from the synthesized movements (see Table 1). Moreover, the highest identification probability from this excerpt was 0.43, which means that it did not clearly identify any other signer from the anonymized movements.

Discussion

This paper shows that simple statistics of the movements of a signer can be manipulated in order to regenerate mocap recordings with a modified identity attribute. The mocap data of SL discourses can undergo various manipulations, such as kinematic identity conversion or anonymization. Moreover, the synthesis algorithm preserves the original temporal structure of the movements, which is crucial because degrading temporal structure could impair the comprehension of the SL discourse.

Up to now, anonymization methods of SL content were modifying appearance, using virtual signers [START_REF] Kipp | Sign language avatars: Animation and comprehensibility[END_REF] or modified videos (e.g., face-swapped videos, where the face of the signer is replaced with another face) [START_REF] Lee | American sign language video anonymization to support online participation of deaf and hard of hearing users[END_REF][START_REF] Bragg | Exploring collection of sign language datasets: Privacy, participation, and model performance[END_REF].

Our technique focuses on controlling the identity in the kinematics of the signers, which could interestingly complement prior approaches in order to provide full anonymity, beyond face or body shape manipulations. Moreover, the proposed algorithm has the advantage that it can render the movements of signers as neutral (i.e., not reflecting the identity of any other signer), by contrast with face-swapping techniques. However, some limitations of the present work should be noted in order to ensure an effective use of these tools in actual applications. First, although we aimed to use SL mocap data as representative as possible of reallife conditions (i.e., spontaneous LSF), the discourses used in the present study were picture descriptions, which may have involved specific linguistic structures more than others (e.g., depicting ones). The different outcomes reported here should be further tested in a wider linguistic context.

Moreover, the present computational findings call for further tests with human participants. Three key problems should be investigated, similarly to prior work on video anonymization [START_REF] Lee | American sign language video anonymization to support online participation of deaf and hard of hearing users[END_REF]: (1) identifiability, by verifying that the ability of human observers to identify the signers is compromised when showing the synthesized modified movements, as compared to the original ones (e.g., with "point-light" displays like in [START_REF] Bigand | Person identification based on sign language motion: Insights from human perception and computational modeling[END_REF] and [START_REF] Troje | Person identification from biological motion: Effects of structural and kinematic cues[END_REF]); (2) comprehensibility, by evaluating the extent to which the observers still understand the SL content in the modified motion examples; and (3) acceptability, by assessing the deaf user perspective on the virtual signers animated with the modified movements and discussing potential use cases (e.g., with focus groups). Should these three fundamental points be validated, the present work could constitute a first step of interest toward automatically controlling the identity of deaf SL users when expressing themselves via virtual signers. Moreover, as shown for videos [START_REF] Bragg | Exploring collection of sign language datasets: Privacy, participation, and model performance[END_REF], preserving anonymity in mocap recordings could increase willingness of SL users to participate in mocap research (e.g., in data collection), which is crucial to develop effective and acceptable technologies.

  , the used markers were (L = left, R = right, F = front, B = back): (1) pelvis, (2) stomach, (3) sternum, (4) LB head, (5) LF head, (6) RB head, (7) RF head, (8) L shoulder, (9) L elbow, (10) LB wrist, (11) LF wrist, (12) LB hand, (13) LF hand, (14) R shoulder, (15) R elbow, (16) RB wrist, (17) RF wrist, (18) RB hand, (19) RF hand. The mocap examples were normalized with respect to size, shape and posture of the signers (see Bigand et al. (2021)). The mocap data of the pelvis marker were ignored as it was set as the origin, which leads to zero vectors. Position and velocity of the body markers were used as temporal features. Velocity was estimated by time differentiation of the mocap position coordinates.

Figure 1 :

 1 Figure 1: The 19 upper-body markers used in the mocap recordings.

Figure 2 :

 2 Figure 2: Distributions of position and velocity data of the RF hand marker along the Z axis, for mocap example 24. Dashed vertical lines represent the means.

Figure 3 :

 3 Figure 3: The two moments of the position and velocity data along the Z axis, for all markers and all 144 mocap examples. Thick lines represent the average statistics of each signer across their 24 examples.

Figure 4 :

 4 Figure 4: The covariance of velocity between body markers (rows and columns) of Signer 2 and Signer 4 in the three dimensions, for mocap example 24. Markers are sorted from the 1 st to the 19 th as presented in Bigand et al. (2021), along X, Y and Z axes. Coefficients correspond to the covariance measures centered and standardized across examples and signers. Blue represent positive covariances, while red represent negative ones. (A) covariance between all markers along the Y axis. (B) covariance between the right hand and arm markers along the Y axis, and the trunk and head markers along the X axis.

Figure 5 :

 5 Figure 5: Schematic representation of the steps used in the synthesis algorithm for the kinematic control of identity.

Figure 6 :

 6 Figure 6: Example of the synthesis results for the identity conversion from Signer 1 to Signer 2 (mocap example 1). Position (up) and velocity (down) data of RF hand marker along the Z axis are shown, for the original mocap recording and synthesized mocap excerpt.

Table 1 :

 1 Output of the automatic signer identification model from original and synthesized movements of Signer 1. The synthesized versions consist of identity conversion into Signer 2 and anonymization. Each output number is the probability that the movements were produced by the signer. Bold numbers represent the highest probability across the six signers.

		Original		Synthesized
			Conversion Anonymization
	Signer 1	0.99	0.00	0.05
	Signer 2	0.00	0.99	0.34
	Signer 3	0.00	0.00	0.14
	Signer 4	0.00	0.00	0.01
	Signer 5	0.00	0.00	0.02
	Signer 6	0.00	0.00	0.43

Acknowledgements

This work has been funded by the Bpifrance (https://www.bpifrance.fr/) investment project "Grands défis du numerique", as part of the ROSETTA project (RObot for Subtitling and intElligent adapTed TranslAtion).


Julien Ancel
email: julien.ancel@polytechnique.edu

Théo Mandonnet
email: theo.mandonnet@ensta-paris.fr

Michel De Lara
email: michel.delara@enpc.fr

French Optimal and Robust Trade Policies Under Carbon Footprint Reduction Constraint and Uncertainties

Keywords: Optimal trade policy, Stochastic optimization, Import shares, Nationally Determined Contribution, Carbon footprint, France

published or not. The documents may come

Introduction

The goal of the first French "Stratégie Nationale Bas Carbone" (SNBC1) was to decrease by 40% territorial greenhouse gases between 1990 and 2030. This strategy did not take into account the carbon footprint from imported emissions, that is, emissions that take place abroad in order to product goods or services consumed in France. Carbon footprint equates indeed to the sum of non-exported territorial emissions and imported emissions. France seeks now to control its carbon footprint, and not only its territorial emissions.

In this context, the French government turned to the "Haut Conseil pour le Climat" (HCC) in 2020 with the following question: "What is the carbon footprint of goods that we import and how to reduce it efficiently?". In its report [START_REF] Le Quéré | Maîtriser l'empreinte carbone de la France[END_REF]), the HCC highlights that imported emissions amounted up to half of the French carbon footprint in 2015. This observation is supported by evidence provided by an OFCE report [START_REF] Maillet | La contribution des émissions importées à l'empreinte carbone de la france[END_REF]). Reducing carbon footprint seems possible, among other solutions, by controlling trade partners; fostering imports from economies with low carbon intensities is a possible lever.

In this paper, we transpose the objectives of SNBC1 in reduction of carbon footprint, instead of territorial emissions, and seek to determine an optimal trade policy for France, under constraints and uncertainty, be they ecological or economic. Constraints emerge from realism purposes, notably from the satisfaction of the French demand of imported goods and services, and from the carbon footprint reduction goal. Uncertainty arises from the trajectories of decarbonization of economies, i.e. from the fulfilling of their Nationally Determined Contributions (NDCs). Taking this uncertainty into account is crucial when it comes to add robustness to an optimal trade policy. This robustness, when compared to optimal perfect foresight trade policies obtained through deterministic optimization, is the key result of our study.

To our knowledge, there is no literature considering optimal trade policy regarding environmental criteria and using import shares as (direct or indirect) decision variables. Most literature dealing with optimal trade policy seems to fit into the context of specific aspects of neoclassical economics, e.g. oligopoly, monopolistic competition and heterogeneous firms, or comparative advantage (e.g. [START_REF] Haaland | Optimal trade policy with monopolistic competition and heterogeneous firms[END_REF], [START_REF] Eaton | Optimal Trade and Industrial Policy Under Oligopoly[END_REF], [START_REF] Jagdish | Optimal trade policy and compensation under endogenous uncertainty: The phenomenon of market disruption[END_REF], [START_REF] Costinot | Micro to Macro: Optimal Trade Policy With Firm Heterogeneity[END_REF], or [START_REF] Costinot | Comparative Advantage and Optimal Trade Policy[END_REF]). Otherwise, most optimization problems under carbon footprint constraints deal with technical systems, especially energy production ones (e.g. [START_REF] Łukasz | Optimal energy planning models with carbon footprint constraints[END_REF]).

The originality of our approach stems from the use of import shares as decision variables in a problem of carbon footprint reduction, and from our explicit stochastic optimization framework.

We now position our contribution with respect to the literature. Diagnoses of the carbon footprints of European countries have underlined the major part of indirect emissions. Similar studies have been done for Germany (Destatis [2019]), the UK [START_REF] Wiedmann | Development of an embedded carbon emissions indicator: Producing a time series of input-output tables and embedded carbon dioxide emissions for the UK by using a MRIO data optimisation system[END_REF]) and France [START_REF] Maillet | La contribution des émissions importées à l'empreinte carbone de la france[END_REF], Le Quéré et al. [2020]). They all emphasize the necessity of taking a consumer point of view to better account for the carbon footprint of consumed goods and in the determination of suitable public policies. Among possible solutions, the control of the carbon intensity of imported goods is highlighted for France in Le [START_REF] Le Quéré | Maîtriser l'empreinte carbone de la France[END_REF]. This approach is supported by previous research [START_REF] Xu | A structural decomposition analysis of the emissions embodied in trade[END_REF]) showing the growing role of imported emissions in carbon footprint of developed countries due to the structure of the value chain. It suggests to better keep track of the carbon impact of imports, and to use trade as a mitigation lever in public policies. [START_REF] Moran | Quantifying the potential for consumer-oriented policy to reduce European and foreign carbon emissions[END_REF] looked at changes in the composition of imports as a consequence of carbon-friendly consumer spending choices. But they only considered individual consumer actions (though massively adopted), whereas our approach takes the point of view of a regulator aggregating at a national level a goal of explicitly controlling imported emissions. Thus, a certain interest for imports choices as a mitigation lever of carbon footprint has been shown both in the literature and by governments, but we did not find previous works effectively determining a trade policy based on a carbon target. In order to determine such policy, we write the optimization problem of a regulator subject to a carbon constraint and explicitly taking shares of French imports represented by each world region as decision variables.

Since we are interested in establishing a long term trade policy, our optimization problem has to take into account not only present but also future carbon intensities of sourcing regions for French imports. These future intensities are widely considered to be uncertain, though countries have provided NDCs after the 21 st Conference Of the Parties (COP) indicating a pledge for their future emissions. den Elzen et al.

[2019] pointed out the uncertainty of major countries (and greenhouse gases emitters) meeting their NDCs. [START_REF] Benveniste | Impacts of nationally determined contributions on 2030 global greenhouse gas emissions: uncertainty analysis and distribution of emissions[END_REF] proposed a quantification of this uncertainty for every country, and assembled estimates of Gross Domestic Products (GDP) and emissions in 2030 and 2050 for all countries produced by the Organisation for Economic Co-operation and Development (OECD), the International Institute for Applied Systems Analysis (IIASA), the Potsdam Institute for Climate Impact Research (PIK) and the "Centre d'Études Prospectives et d'Informations Internationales" (CEPII). Thus we choose an optimization framework where this uncertainty appears explicitly. We will see that this approach leads to more robust trade policies than policies obtained under a perfect foresight approach, where the future emissions and GDPs are supposed known in advance.

This paper is organized as follows. In Section 2, we explicitly present our notations, then link them to the available data to finally formulate the stochastic optimization problem corresponding to least trade cost minimization under constraints and uncertainties. In Section 3, we display different numerical results; we compute optimal robust trade policies for 36 trading partners, also aggregated in 3 and 2 trading world regions; we discuss the interest of our stochastic approach. We conclude in Section 4.

Data and method

The problem studied here is the determination of an optimal French trade policy in order to follow the recommendation of the HCC in terms of carbon footprint reduction. This goal is a 40% decrease between 1990 and 2030. The framework used in this study is the so-called two-stage stochastic optimization [START_REF] Shapiro | Lectures on stochastic programming: modeling and theory[END_REF]). The two stages correspond to the years 2015 and 2030. Decisions taken in the first stage (2015) are the import shares from each region. On the second stage, if the footprint target is not reached, France has to buy emissions quotas in order to cover the difference between its footprint and the target.

The random aspect of the problem is to be found in the French demand in 2030, in the carbon intensities of world economies in 2030 and in the price of the European carbon emission quotas in 2030.

The monetary amount of quotas France has to buy to cover an excess of carbon footprint over the reduction target is our penalization criterion. The estimation of the distribution of the random variables presented in the following of this section is based on the work of [START_REF] Benveniste | Impacts of nationally determined contributions on 2030 global greenhouse gas emissions: uncertainty analysis and distribution of emissions[END_REF]. This study generates GDP and emission data for all countries in 2030 based on Nationally Determined Contributions announced by the countries on the occasion of Paris COP 21 and revised for COP 26.

GDP and greenhouse gases emissions in 2030

Data used in our study comes from the work of [START_REF] Benveniste | Impacts of nationally determined contributions on 2030 global greenhouse gas emissions: uncertainty analysis and distribution of emissions[END_REF]. In the latter paper, the authors gather data and estimations of GDP and greenhouse gases (in CO 2 eq) for every country, for the years 2015, 2030 and 2050. These data and estimations are furnished by 4 organisms: CEPII, IIASA, OECD and PIK. For future years, estimations of the distribution of emissions based on NDCs are proposed by each of these organisms, based on Nationally Determined Contributions given by the states themselves. When no NDC is given by a country, the organisms estimate themselves, as experts, the evolution of emissions of this country.

As a matter of fact, these estimations are done 5 times. Indeed, along with the IPCC reports it is commonly accepted, in the literature dealing with climate change and its attenuation, to consider 5 stylized (and standardized) possible versions of the world and its evolution: these are the Shared Socioeconomic Pathways (SSP). These SSPs account for the diversity of consequences of climate change and of climate policies. They are used by [START_REF] Benveniste | Impacts of nationally determined contributions on 2030 global greenhouse gas emissions: uncertainty analysis and distribution of emissions[END_REF] to build trajectories of greenhouse gases emissions associated with different sets of climate policies. 5 SSPs are considered: Sustainability (SSP1), Middle of the Road (SSP2), Regional Rivalries (SSP3), Inequalities (SSP4) and Fossil-fueled Development (SSP5). Uncertainty arises from NDCs, for example on the capacity of a country to meet its own NDC or on the reliability of the design of the NDC.

Thus [START_REF] Benveniste | Impacts of nationally determined contributions on 2030 global greenhouse gas emissions: uncertainty analysis and distribution of emissions[END_REF] build distributions of emissions for each country for each SSP and for each GDP data source. In our study, we take advantage of this work and build 20 scenarios of carbon intensities of world regions. A given scenario of ours is characterized by the choice of a tuple (SSP, GDP data source) -e.g. (SSP1, CEPII). For every region, we compute the mean of the distribution of emissions produced by [START_REF] Benveniste | Impacts of nationally determined contributions on 2030 global greenhouse gas emissions: uncertainty analysis and distribution of emissions[END_REF] for this region for this tuple (SSP, GDP data source). Finally, we divide this mean by the GDP of the region given by the selected GDP data source. This yields a carbon intensity of each region for this tuple (SSP, GDP data source).

Here below, we detail the mathematical formalization.

Notations for the decision problem

In order to write the two-stage optimization problem of a regulator -seeking to reach a carbon footprint target by modifying French trade policy, while minimizing the value of emissions quotas bought -we take the following notations, in which stage t = 0 refers to year 2015 and stage t = 1 refers to 2030.

As seen in §2.1, we use the estimations gathered by [START_REF] Benveniste | Impacts of nationally determined contributions on 2030 global greenhouse gas emissions: uncertainty analysis and distribution of emissions[END_REF] to derive 20 scenarios (one per couple SSP×organism furnishing data) in order to account for uncertainty arising from NDCs. We envision this set of scenarios to explore several possibilities of evolution of the variables presented below. More generally, scenarios are supposed to form a finite set S, with the generic scenario being referred to as s, and having probability π s ≥ 0 (with s∈S π s = 1). In the equiprobable case, the probability of any scenario s is π s = 1/|S|, where |S| is the cardinality of the finite set S of scenarios.

The world regions with which France trades constitute the finite set R. The generic trading partner is referred to as r. We consider a set of 36 sourcing regions that gather 90% of the French imports in 2015, according to the French "Direction Générale des Douanes"1.

Let I 0 be the vector of French demand fractions, that are imported by France from each world region r ∈ R at stage t = 0, constituting the first decision variable as

I 0 = {I r 0 } r∈R . (1)
These decisions are taken at stage t = 0 (2015), hence the lower index in I 0 . Note that the region "France" appears in the vector I 0 in (1), so that the sum of the I 0 coefficients equates to 1. The fraction associated with the region France is blocked at its 2015 level: we do not account for the possibility of relocation, nor of offshoring. Let C s 1 denote the vector of carbon intensities of the different world regions r ∈ R at stage t = 1 (2030), in the scenario s:

C s 1 = {C r,s 1 } r∈R . (2)
1https://lekiosque.finances.gouv.fr/site_fr/telechargement/telechargement_conjoncture.asp

These uncertainties are supposed to be revealed at stage t = 1, hence the lower index in C s 1 . The different scenarios s ∈ S allow to account for uncertainty in NDCs implementations.

Let E 1 = 249 MtCO 2 eq denote the 2030 carbon footprint target, which corresponds to a decrease of 40% from the French carbon footprint of 1990. Let D s 1 be the total French demand of importations at stage t = 1, which depends on the scenario s. Finally, we denote by Q s 1 the quantity of emissions quotas bought by France at stage t = 1 and by P s 1 its unitary price at stage t = 1, in the scenario s. France has to buy emissions quotas if its target E 1 is not reached.

From data to ingredients for the decision problem

Estimates of the French GDP yield D s 1 , the value of the French demand of importations at stage t = 1 in scenario s. Imports are considered to represent the same part of D s 1 at t = 1 as at t = 0, that is, 30% according to INSEE. This is a strong hypothesis since this part has varied significantly in the past decades, but it enables to isolate the effects of the modification of the trade policy on carbon footprint. Denoting by gdp F R,s 1 the French GDP in 2030 in scenario s, we thus set

D s 1 = 0.3 × gdp F R,s 1 . (3)
On the other hand, intensities C r,s 1 are obtained by dividing, for each region r, the GDP by the greenhouse gases emissions. Denoting by ghg r,s 1 the emissions of region r in 2030 in scenario s, C r,s 1 is calculated as

C r,s 1 = ghg r,s 1 gdp r,s 1 , ∀(r, s) ∈ R × S . (4)
The price of a carbon quota is taken according to the conclusion of Commission Quinet [START_REF] Quinet | La valeur de l'action pour le climat -Une valeur tutélaire du carbone pour évaluer les investissements et les politiques publiques[END_REF]) which determined a long term carbon price scenario to use in France. Thus, the price of carbon is fixed at 250€ 2018 per ton in 2030. This price is also in line with the approximate levels reached by the price of an emission allowance on the EU-ETS permits market on the long run due to cap restriction.

To strengthen the economic realism of the derived trade policies, we suppose that, for each world region r ∈ R, the demand fraction I r 0 is bounded above by

I r 0 = δ × I r,real 0 , ∀r ∈ R , (5)
where I r,real 0 is the real French import share from region r in 2015, and where the value of δ is set to 1.6 as explained below. The term δ limits the accessible variation in commercial exchanges between France and a region, by encompassing all possible costs of such variation that are not explicitly modeled in this study. These costs may arise from various reasons such as historical or geopolitical links between the region and France, geographical distance or transport costs. It also reflects the limit of each region in production of goods and services that it may export: France cannot possibly import more from a region than what it offers to export. Finally, the parameter δ may also be seen as one of relative sensitivity to the carbon constraint: an increased parameter increases the relative weight of the respect of the carbon constraint among the factors that determine French trade policies. The arbitrary value of δ = 1.6 is the result of several tries. Lower values (< 1.4) contract indeed the set of admissible policies very close to the current policy. These values are relevant in worlds that do not care about carbon constraints. Since this study seeks precisely to measure the impact of caring about imported carbon content when trading, low values of δ are not suitable. Much higher δ values (> 1.8) yield trade policies where France exchanges with countries that could not realistically export that much goods, or, in some scenarios, where France leaves completely its historical partners. These high values make the carbon constraint prevail on other determinants of trade and, more importantly, on the production limits of sourcing regions. Thus, higher values are not relevant either, and we keep δ = 1.6 which balances between a relative importance of the carbon constraint in the determination of trade policies, and realistic variations from the 2015 situation with a 15 years horizon.

Formulation of a stochastic optimization problem

We consider 20 scenarios, one per couple (SSP × organism furnishing data), all equiprobables. With the notation in §2.2, this means that |S| = 20 and that the probability of scenario s is π s = 1/|S| = 1/20.

First, we write the problem formulation in an extensive linear form, with recourse variables {Q s 1 } s∈S in addition to the original decision variables I 0 in (1). Second, we write an equivalent problem in the original decision variables I 0 , having the property to be a convex program.

Extensive linear form (in initial and recourse decision variables)

With the notations in §2.2, we write the optimization problem of a regulator as

min {I r 0 } r∈R ,{Q s 1 } s∈S s∈S π s P s 1 Q s 1 (6a) 0 ≤ Q s 1 , ∀s ∈ S , (6b)
0 ≤ I r 0 ≤ I r 0 , ∀r ∈ R , ( 6c
) r∈R I r 0 = 1 , (6d)
D s 1 r∈R I r 0 C r,s 1 -Q s 1 ≤ E 1 , ∀s ∈ S . (6e)
In the optimization problem ( 6), the regulator seeks to minimize the expected cost (6a) of buying emissions quotas at stage t = 1 (2030), under constraints that we comment right now. Constraint (6b) is a constraint of nonnegativity of the quantity Q s 1 of emissions quotas. Constraint (6c) corresponds to the nonnegativity of the French demand fractions I r 0 and to a maximal variation of I r 0 from the 2015 situation (t = 0). This constraint of realism accounts for all the links between countries not explicitly taken into account in our study that also influence trade between France and these countries. This maximum is a fixed parameter detailed in (5): each region can export to France up to 1.6 times the part of French demand it represents at stage t = 0. Constraint (6d) accounts for the satisfaction of French demand. Constraint (6e) represents the achievement of the carbon footprint target E 1 for all scenarios s ∈ S, thanks to the possibility of recourse by means of additional emissions quotas Q s 1 . Finally, note the difference in available information when the two decision variables are determined. I 0 is deterministic; this reflects a decision taken at stage t = 0. On the contrary, the decision Q s 1 depends on the scenario s because it is a function of French demand and carbon intensities of regions at t = 1, which all depend on the scenario s.

Convex form (in initial decision variables)

In case the price P s 1 of emission quotas bought by France at stage t = 1 is nonnegative for any scenario, that is, when

P s 1 ≥ 0 , ∀s ∈ S , (7)
another interesting formulation is possible. Using the form of the objective function in (6a) and the constraints (6e), it is easy to observe that the variable Q s 1 is, at the optimum in the optimization problem (6), fully determined as a function D s

1 r∈R I r 0 C r,s 1 -E 1 +
of import shares, future French demand and future carbon intensities, where x + = max{x, 0}. Thus, in this setting, the optimization problem (6) rewrites as

min {I r 0 } r∈R s∈S π s P s 1 D s 1 r∈R I r 0 C r,s 1 -E 1 + (8a) 0 ≤ I r 0 ≤ I r 0 , ∀r ∈ R , ( 8b
) r∈R I r 0 = 1 . ( 8c
)
In this form, the objective function (8a) of problem ( 8) is clearly a convex function of the vector of import shares {I r 0 } r∈R thanks to the polyhedral convex (• • •) + terms in (8a). Constraints (8b)-( 8c) are linear and depend also only on variables {I r 0 } r∈R . Thus, this form places problem (8) in the framework of convex optimization, with the property that solutions are then expected to be (generally) inside their constraints domain (whereas solutions of linear programs lie on the border). Thanks to such convex formulation, we obtain robust solutions for the import shares {I r 0 } r∈R .

Numerical results

Here, we solve numerically the optimization problem (8). For this purpose, we keep, for each country, the mean of its emissions and GDP in 2030 given by the 4 organisms (OECD, IIASA, PIK, and CEPII) for each SSP. This yields 20 equiprobable scenarios (4 organisms × 5 SSPs), that are input data for the convex optimization problem (8).

Once obtained the optimal solution {I r 0 } r∈R -that we coin "robust optimal policy" (more evocative than stochastic optimal policy) or "robust (stochastic) optimal policy" -we compare its performances with other pre-determined policies:

• the current trade policy;

• the 20 "perfect foresight policies" that are solutions of problem ( 6) when all the probabilities π s , in the family {π s } s∈S , are set to zero except for one scenario.

Let us be more specific regarding this last case. Condider a single scenario s ∈ S (for instance, scenario SSP1 estimated by IIASA). The perfect foresight policy associated with scenario s solves the following deterministic program:

min {I r 0 } r∈R ,Q s 1 P s 1 Q s 1 (9a) 0 ≤ Q s 1 , ( 9b
)
0 ≤ I r 0 ≤ I r 0 , ∀r ∈ R , ( 9c
) r∈R I r 0 = 1 , (9d)
D s 1 r∈R I r 0 C r,s 1 -Q s 1 ≤ E 1 . ( 9e
)
This is problem (6) with π s = 1 (e.g. π IIASA-SSP1 = 1) and all other π s = 0, that is, we optimize thinking that the future will be scenario s for sure. These perfect foresight policies are solution of deterministic optimization problems, as the scenario is known in advance at stage t = 0, that are linear programs; as a consequence, they yield "corners" solutions, that are less robust to hazard. Indeed, resulting import shares are closely associated with the scenario they are optimized for. For pessimistic scenarios, this gives policies concentrated in few countries, thus closely associating the success of this policy -in reaching the carbon footprint target -to the successful implementation of NDCs of these few countries. Other scenarios may be very optimistic as they feature the completion of NDCs for countries that, in most other scenarios, fail to reduce their emissions. Optimistic scenarios yield optimal policies that rely, at least partially, on such countries. Thus resulting policies may also be considerably stretched between all world regions. In a policy approach, the latter optimistic case is more detrimental. Indeed, if France imports from fewer countries than in the current situation, and estimates that, in the end, these partners might not sufficiently reduce their emissions -leading to supplementary needs for emissions quotas -France may decide to help them, if the (social) cost of helping is lower than the final cost of supplementary emissions quotas and climate related damages. On the contrary, such recourse is not an option when the decided trade policy has substantially diversified regions of imports. This policy trap -bad trade policy choice and no possible recourse -exists because of the cherry picking of the optimal solution when deterministic optimization along a particular scenario is done. The robust (stochastic) optimal policy, that we advocate for, avoids this issue by taking all scenarios into account altogether.

Results for 2 aggregate regions

To better perceive the robustness of decarbonization policies obtained through stochastic optimization, we present a 2D case in Figure 1, where only 2 regions, apart from France, are considered: Europe and the rest of the world. It can easily be seen that all the deterministic optimal solutions saturate some constraints -they are all near corners -but the robust (stochastic) optimal solution is not in a corner of the admissible polygon (a segment, here), which makes it more robust to hazard. The current situation is clearly distinct from any of the presented optimal policies (deterministic or stochastic), but not that far from the robust (stochastic) optimal policy. In the robust (stochastic) optimal solution, the European import share is increased from 70% to roughly 80% compared to the current situation. This reflects that, even when uncertainty is accounted for, European economies are significantly more decarbonized than those from the rest of the world. However, this general trend does not apply to all scenarios, and this is why there exist deterministic optimal solutions with a decreased European share when compared to the current situation: these solutions are optimal in scenarios where the trend is not verified. In any case, the robust (stochastic) optimal policy leads to a more balanced, more robust trade policy.

Results for 3 aggregate regions

Going on, we simplify the original problem with 36 countries to only 3 big world regions: Europe, America and Asia. We then represent, in R 3 , the import shares for the robust (stochastic) optimal policy, the current policy and the 20 perfect foresight policies. In Figure 2, we can see that the 20 "deterministic optimal solutions" (perfect foresight) are spread wider in [0, 1] 3 because they lie on corners of the polygon of feasible solutions, as optimal solutions of a linear program. This explains the robustness of the robust (stochastic) optimal policy: the spread of the deterministic optima reveals the diversity of possible futures. Note that, as explained above, some perfect foresight solutions increase the balance between import shares for all 3 regions, even though it appears that, in most scenarios (namely all scenarios except the few these optima are optimized for), Asian and American emissions tend to be significantly higher than European ones. Thus, the "wisest" option may not necessarily be to balance between regions. We insist therefore on the fact that perfectly adapting to only one future may be detrimental if another future is in fact realized, whereas the robust (stochastic) optimal policy balances all of these possible futures.

Results for 36 world regions

Results with 36 world regions illustrate how the robust (stochastic) optimal policy is more robust to randomness than other pre-determined policies. On the one hand, the optimal policies for a deterministic problem fare better than the optimal solution of the stochastic problem for the particular futures they were designed for. On the other hand, they are also less robust to randomness: their frequencies of necessary recourse to additional carbon quotasi.e. of failure to meet the carbon footprint target (event Q 1 > 0)are significantly greater than what is obtained with the optimal solution of the stochastic problem (Figure 3). In every scenario, the ratio between necessary emissions permits and imported carbon footprint amounts at least to 50%, whatever the policy is. However, this ratio is systematically the lowest for the robust (stochastic) optimal policy. Thus, this policy leads not only to the lowest probability of having to buy emissions permits (Figure 3) and the most on the left distribution of final carbon footprint , but also to the lowest relative importance of emission permits in the meeting of the carbon footprint target.

Figure 3: Comparison of frequencies of the event "No necessary supplementary quota" (i.e. Q 1 = 0) for the robust (stochastic) optimal policy, the current and mean policy and the 20 perfect foresight policies We look at the actual import shares in the 36 regions framework (Figure 4), as they currently amount to 90% of French imports (French "Direction Générale des Douanes" data). The robust (stochastic) optimal solution is characterised by higher shares devoted to European countries, notably Germany, Spain and the UK. It is also noteworthy that the robust optimum policy leads to zero imports from China, Belgium, or Poland, for instance. This notable difference with the current policy is explained by the projections of the carbon intensities of these countries (see Figure 5 in Appendix A). Indeed, the robust (stochastic) optimal policy concentrates on the countries that have the lowest carbon intensities across the different SSPs. This entails the disappearance of China and Poland -that have clearly higher carbon intensities -while we observe a significant increase of import shares from Germany, Italy, the United States, Spain, Netherlands, the United Kingdom, or even Switzerland. All those latter countries have the relatively lowest carbon intensities under uncertainty. In particular, Japan and Switzerland do not display GHG emissions uncertainty, as experts are very confident in these countries implementing successfully their NDCs (see [START_REF] Benveniste | Impacts of nationally determined contributions on 2030 global greenhouse gas emissions: uncertainty analysis and distribution of emissions[END_REF] for an explanation). Note that Belgium, that has a slightly higher carbon intensity distribution, also drops almost to zero like Poland and China. As emissions distributions of countries may clearly be separated in different classes, the robust (stochastic) optimal policy concentrates imports in regions belonging to the class of low emissions countries, as these countries are more carbon-efficient in all considered futures.

Conclusion and discussion

Using import shares as decision variables -and aiming to minimize the volume of carbon quotas bought by France in 2030 -we have seen, in this study, that an explicitly stochastic optimization approach would lead to a more robust trade policy than any deterministic one. When tested in various scenarios, the former policy is less likely than the latter to miss the carbon footprint target -i.e. to necessitate the buying of emissions permits -and, when it does, fewer permits are required. On the contrary, perfect foresight (that is, anticipative) policies built to be optimal in a certain version of the future behave better in this version than the robust solution, but worse in all other futures. The greater robustness of our stochastic optimal policy advocates for the diffusion of stochastic optimization frameworks when dealing with scenarios of future, rather than optimizing according to each scenario.

It appears, in Figure 4, that the current French policy is relatively close to the robust (stochastic) optimal policy. This may be explained by the fact that France imports essentially from the European Union -which is composed of countries that have, on the whole, lower carbon intensities of production than the rest of the world. The distributions of GHG emissions provided for these countries by IIASA, CEPII, OECD and PIK have also generally lower spread. It reflects the relative confidence of the experts in, at least, the design of their NDCs and, perhaps, in their strategies to fulfill their NDCs. These lower centers and lower spreads of distributions may explain why countries (mostly non-European) that represented low shares of French total imports (namely less than 1%) or countries subject to debate between experts (such as Belgium) disappear in the robust (stochastic) optimal policy (Figure 4).

It is noteworthy that trade policies are often determined through border tariffs, and not directly import shares. The method we propose could be coupled with a general equilibrium calculation to derive the optimal border tariffs associated with these import shares. Perhaps the trade policy could also be influenced by non-tariff measures, such as suggested in [START_REF] Ederington | Nontariff Measures and the World Trading System[END_REF]. Finally, the economic realism of this study is obviously debatable: a real trade policy has to respect WTO agreements and, more generally, to fit in a complex geopolitical and economic context. Still, the objective is to encourage the use of explicitly stochastic frameworks to make decision-making more robust, rather than to yield an effective French optimal trade policy.

Figure 1 :

 1 Figure 1: Scatter plot of robust (stochastic) optimal policy, current policy and the 20 perfect foresight policies (2 aggregate regions)

Figure 2 :

 2 Figure 2: Scatter plot of robust (stochastic) optimal policy, current policy and the 20 perfect foresight policies (3 aggregate regions)

Figure 4 :

 4 Figure 4: Import shares from the current 36 biggest suppliers of France in the current policy (bottom left in blue), robust (stochastic) optimal policy (bottom right in red) and a selection of 2 (among 20) perfect foresight policies (upper, in black)