Yun Lee Fabian

Shaiful Lee

Hashim Jahari

Samsudin Khairulmizam

Fakhrul Noor 'ain Kamsani

Rokhani Zaman

Wei Yang

Lim

Earthquake Prediction using Time-Frequency Distribution Methods

Introduction

Earthquakes are generally a product of seismic movements, namely by tectonic plates abrupt movements, releasing waves of energy. These waves that travel through the crust of the earth are the cause of the shaking that humans feel. For earthquakes with higher magnitudes, civil structures may be displaced laterally or even collapse [START_REF] Tao | Time-Frequency Energy Distribution of Ground Motion and Its Effect on the Dynamic Response of Nonlinear Structures[END_REF]. For countries that are more earthquake-sensitive, experience has driven them to innovation for countermeasures, typically the development of vibration control devices to increase the resilience of civil structures [START_REF] Li | On vibration suppression and energy dissipation using tuned mass particle damper[END_REF][START_REF] Spencer | State of the Art of Structural Control[END_REF].

An early and accurate prediction of earthquake occurrences may lead to better disaster management and preparations for structural damage control, and ultimately, leads to reducing injuries and deaths caused by earthquakes, as well as economic losses. As such, earthquake prediction is thought to be the 'holy grail' of seismology, as it is near impossible given the various factors and margin of error in time, place, and magnitude [START_REF] Hough | The Great Quake Debate: The Crusader, the Skeptic, and the Rise of Modern Seismology[END_REF].

As the spark for this research, Turkey and Syria saw an earthquake doublet: magnitude 7.8 occurring on 6th of February 2023, saw a total casualty of over 52,000 deaths in both Turkey and Syria [START_REF] Naddaf | Turkey-Syria earthquake: what scientists know[END_REF][START_REF] Dal Zilio | Earthquake doublet in Turkey and Syria[END_REF], and approximately 9 hours later a magnitude 7.6 earthquake occurred north-northeast from the former [START_REF] Dal Zilio | Earthquake doublet in Turkey and Syria[END_REF]. The estimated total structural damages and immediate production capacity losses are upwards of US$54 billion [START_REF] Demiralp | The Economic Impact of the Turkish Earthquakes and Policy Options[END_REF].

As these seismic movements or waves can be measured and recorded, with the intensity of the waves noted as magnitude, the recorded waves are essentially a time-series function describing the occurrences of earthquakes and their respective magnitudes. Time-Frequency Distributions (TFD) are used in signal processing to analyse and characterise signals in the time and frequency domains simultaneously, and have been used to analyse earthquake motion measured using accelerograms [START_REF] Demetriu | Time-frequency representations of earthquake motion records[END_REF][START_REF] Huerta-Lopez | Time-frequency analysis of earthquake records[END_REF]. In this research, we propose to analyse earthquake records using TFDs method to find precursory patterns in the time-frequency domain for earthquake prediction.

Methodology Data Sources

Datasets of the Turkish (2023), Indonesian (2018), Filipino (2022), and New Zealander (2016)r earthquakes were used for the analysis conducted in this research [START_REF] Dincer | EarthQuake Turkey 06 Feb 2023 -Latest[END_REF][11][12][13]. For this research, we will only be focusing on the timestamp and magnitude data provided in the dataset, as the time-series function to undergo Time-Frequency Distributions (TFDs). 1), all datasets appear to show irregular sampling, which is not meeting a crucial criteria for applying the key function for TFDs -Fast Fourier Transform (FFT) -i.e. the time-series data provided must be, or is assumed to be, regularly sampled [START_REF] Cochran | What is the fast Fourier transform[END_REF]. Therefore, a data interpolation method is used, to estimate and produce a regularly sampled data based on each dataset. The interpolation method used is linear interpolation, which is done by estimating values between two adjacent points using a straight line. As earthquake occurrences are not a high frequency phenomena, losing high frequency information due to linear data interpolation can be considered minimal impact to analysis result reliability, if not negligible.

Data Interpolation

Time-Frequency Distributions (TFDs)

Time-frequency analysis is mainly achieved by formulating a function which represents the input signal in both time and frequency spectrum simultaneously, instead of either within time or frequency domain. This function or process can then be changed to suit the various use-cases that contain a time-varying spectrum [START_REF] Hlawatsch | Time-Frequency Analysis[END_REF]. There is no objective representation for signal energy changes over time and frequency, and this resulted in several proposals [START_REF] Hlawatsch | Time-Frequency Analysis[END_REF].

Short-Time Fourier Transform (STFT)

The Short-Time Fourier Transform (STFT), also known as spectrogram distribution, is a signal processing method used to analyse multi-component non-stationary signals.

The key idea of this method is to apply FFT over a fixed width window moving over time, resulting in the change of spectral information of a signal over each frame.

Wigner-Ville Distribution (WVD)

The Wigner-Ville distribution (WVD) for a continuous signal is given as follows:

𝑊𝑉 𝑡,ω = 1 2π ∫ 𝑠 𝑡 + τ 2 ( ) 𝑠 * 𝑡 -τ 2 ( ) 𝑒 -𝑗ωτ 𝑑τ (1)
With its kernel function being 1.

The WVD covers most mathematical properties involved in analysis; real-valued, maintains time and shifting frequencies, and satisfies marginal features that generate an instantaneous power and power spectrum.

Choi-Williams Distribution (CWD)

One of the most significant and well-known representations of signal energy distribution is the quadratic class, or Cohen's class [START_REF] Li | Analysis of time-frequency energy for environmental vibration induced by metro[END_REF]. However, Cohen's class distribution is hampered by unwanted cross-terms between frequency components [START_REF] Pereira | Reducing cross terms effects in the Choi-Williams transform of mioelectric signals[END_REF]. Among the best methods for reducing these terms to the minimum is the Choi-Williams [START_REF] Choi | Improved time-frequency representation of multicomponent signals using exponential kernels[END_REF] class kernel (CWD), with the discrete form expressed as follows:

𝐶𝑊𝐷 𝑛,ω = 2 𝑘= 𝑁/2 𝑁/2-1 ∑ 𝑒 -𝑖2ω𝑘 𝑚=-𝑀/2 𝑀/2 ∑ σ 4π𝑘 2 𝑒 - σ 𝑚-𝑛 ( ) 2 4𝑘 2 𝑥 𝑚 + 𝑘 ( ) 𝑥 * 𝑚 -𝑘 ( ) (2)
Pereira proposed more cross-term reduction for Cohen's class so that only the positive matrix part contains the data information, while the negative part is related to the cross-terms and therefore excluded.

Mean Energy Concentration

The mean energy concentration of TFDs refers to the distribution of energy across the time and/or frequency axes in the TFD representation of a signal, and provides insights into the dominant time/frequency/time-frequency regions where the signal's energy is concentrated. For this research, we will be focusing on the mean energy concentration on the time axis of the TFD representation of the earthquake datasets, to observe the dominant time regions of the earthquake datasets for any precursory behaviour of major earthquakes in the time domain. The mean energy concentration of the discrete time-frequency is calculated, via determining time slots of the TFD with a marginal magnitude (N). For time slot N, and number of samples S:

𝑀 𝑡,ω = 1 𝑁 𝑡=𝑁 𝑆 ∑ 𝑓 ∑ 𝑇𝐹(𝑡, ω) (3)
The result from (3) is then normalised as follows:

𝑀 𝑛𝑜𝑟 = 𝑀 𝑡,ω -𝑀𝑖𝑛[𝑀 𝑡,ω ] 𝑀𝑎𝑥[𝑀 𝑡,ω ] -𝑀𝑖𝑛[𝑀 𝑡,ω ] (4)

Results

The TFDs chosen for this research will be compared in regards to their accuracy and energy concentration as they are utilised to represent and analyse the earthquake data.

The earthquake dataset of Turkey was obtained from Kaggle, which spans from 4th of January 2023 to 6th of February 2023 (the onset of the reported major earthquake).

The data is then linearly interpolated to preserve the low frequency (sub-Hz) content, to produce a dataset with constant sampling rate of 200 seconds/sample. This allows the interpolated data to be analysed using the time-frequency analysis methods STFT, CWD, and WWD, as a constant sampling rate of a signal is a prerequisite for any frequency or time-frequency analysis.

The sampling rate of 200 seconds is chosen from the average of the sampling difference in the original dataset. The normalised mean energy concentration is calculated using a time slot of 10.

The interpolated data is then compared to each set of heatmap plots and the normalised mean

Discussions

From the TFD visualisations of the earthquakes, it can be said that the TFDs were able to pinpoint and provide more insight on the onsets/primary/major of the earthquake incidents, providing the spectral information of the earthquakes, as well as the frequency variations over time.

On observation of the comparisons of the earthquake data, heatmap of respective TFDs, and the mean energy concentration derived from the respective TFDs, cropped to a window containing the major earthquakes (Figures 1a-3a for Turkish Earthquake, Figures 1b-3b for the Indonesian Earthquake, Figures 1c-3c for the Filipino Earthquake, and Figures for the New

Zealander earthquake 1d-3d), it is found that there are a few minor peaks visualised in the mean energy concentration graphs, which were initially barely, if not unnoticeable, in the earthquake signals, raw or interpolated.

The resulting elapsed time in seconds of the peaks since the first sample timestamp of the respective datasets were noted, and the estimated timestamps in UTC are calculated. It is noted that the estimated timestamps are considered estimated, as the mean energy concentration graph were calculated in time bins of 5 or 10 interpolated samples depending on the dataset, resulting in a mean energy concentration value representing windows of 2000 seconds (Turkey, New Zealand), 13000 seconds (Indonesia), and 7500 seconds (Philippines). On comparison of these estimated timestamps with the original datasets, the nearest local peak magnitudes are noted, and a summary of each minor peak corresponding to minor earthquakes prior to each major earthquake is tabulated in Tables 2345.

In the analysis of the Indonesian Earthquake, it is noted that CWD observed a major peak towards the timestamp equivalent to 29th of September 2018 (see Figure 2b, as compared to Figures 1b and3b). This is the result of endpoint effect which causes distorted or inaccurate representations of the signal's time-frequency content near the edges of the signal segments (29th September being the end of timestamp range) [START_REF] Lv | Optimization of the End Effect of Hilbert-Huang transform (HHT)[END_REF]. According to reports, it is also possible that this major peak stems from the aftershocks of the major earthquake which lasted until 30th

September 2018 [START_REF]Map of epicentres of mainshock and aftershocks in first 24 hours[END_REF]. Regardless, for the purpose of conformity and summary, this peak will be treated as an anomaly, and the mean energy concentration of the Choi-Williams TFD will be re-normalised in the following table ( From the tables documenting the observation of the minor earthquake peaks, it is shown that for all major earthquakes in question, there are at least 1 minor earthquake occurring between 20

to 70 hours prior. The minor earthquakes are noted to be of several magnitudes lower than the major earthquakes. On the raw earthquake data, it is noted that this lower magnitude level (4.0 to 6.0) is somewhat common (Table 6). Typically, magnitudes 4.0 to 5.0 are barely noticeable, and magnitudes 5.0 to 6.0 are noticeable but incidents or injuries caused by these magnitudes are extremely rare [START_REF] Rafferty | Richter scale. Encyclopedia Britannica[END_REF]. Observing the number of minor earthquakes recorded against the number of major earthquakes recorded, the ratio of recorded minor earthquakes (magnitudes 4.0 to 6.0) to major earthquakes (magnitudes 6.0 and above) are staggeringly high (87.9:1, 301:1, 51:1, 73:1 for the Turkish, Indonesian, Filipino, and New Zealander earthquakes respectively). This indicates that the occurrences of minor earthquakes alone does not constitute a precursor to the occurrence of major earthquakes, and as such, observing the time-series data of earthquake signals is insufficient to predict major earthquakes.

Magnitude

In support of this research's proposed method, as shown in Figures 123, the mean energy concentration of the time-frequency distribution data of the earthquake signals provides possible hints towards the major earthquake occurrences, with peaks corresponding to some of the minor earthquakes. However, it is also noted that not all recorded minor earthquakes have a corresponding mean energy concentration peak. In other words, the occurrences of minor earthquakes within a time period is not the sole factor in contribution to the mean energy concentration of that time period, but the mean energy concentration is indicative of a signal or data pattern which could indicate a major earthquake.

Limitations and Further Research

This research employs only the magnitude of earthquakes in the investigated regions and their respective time detected to form a time-series signal in the method of analysis, and therefore, is able to predict only an estimated time, magnitude, and the region of major earthquake (time, size, and location) [START_REF] Wyss | Why is earthquake prediction research not progressing faster?[END_REF]. Various other measured parameters have yet been considered and so, the authors exhort for further research to investigate the potential of various parameters such as GPS coordinates and depth, in contribution to the accuracy of earthquake prediction.

The time window between the precursor pattern to the major earthquake found in this research ranges between 20 to 70 hours, which perhaps is acceptable for a prediction timeframe, providing enough time for institutions to warn the general public, however the authors believe that the accuracy of the time window can be improved.

Conclusion

In this study, we have explored the utilisation of Time-Frequency Distributions (TFDs) on seismic signals for earthquake prediction, demonstrating the effectiveness of TFD-based techniques and enabling the identification of patterns of minor earthquakes leading up to major earthquakes.

Through extensive experimentation on real-world earthquake datasets, the results indicate that TFD-based techniques show great potential in accurately predicting seismic events. The incorporation of TFDs provides valuable insights into the underlying dynamics of earthquakes, while the identification of subtle precursor patterns enables improved prediction capabilities.

This information can aid in the development of advanced prediction models and early warning systems, and consequently, proactive measures for disaster (earthquake) management and potentially reduce human and economic losses caused by earthquakes. As this research merely peeks into the usage of TFDs to analyse earthquake signals, further investigations are needed to explore its full potential in predicting earthquakes. Future studies should consider incorporating additional data sources and parameters, such as geodetic measurement, atmospheric data, latitude, longitude, altitude, and depth measurements, to enhance the predictive capabilities. Moreover, the development of robust and scalable TFD algorithms will facilitate real-time analysis and improve the efficiency of earthquake prediction systems.

In conclusion, the analysis of earthquake data using TFDs opens up new possibilities in earthquake prediction, a feat which many considered to be impossible previously. The advancements made in this study contribute to the growing body of knowledge in seismology, offering valuable insights for researchers and practitioners working towards improving our understanding of earthquakes and developing effective strategies for disaster damage control.

  energy concentration plot of the TFDs (spectrogram/STFT distribution, Choi-Williams distribution, Wigner-Ville distribution).The entire process above is repeated using Indonesia's 28th September 2018 earthquake dataset (dataset spanning from 1st September to 30th September 2018), Philippines' 27th July 2022 earthquake dataset (dataset spanning from 1st July 2022 to 31st July 2022), and New Zealand's 16th November 2016 earthquake dataset (dataset spanning from 1st November 2016 to 30th November 2016). It is noted that the datasets for these earthquakes have varying ranges of samples, and therefore, several parameters are adjusted, namely the resampling rate for interpolation of data, the window sizes for the TFDs in order to optimise both time and frequency resolution for STFT, and the time slots used for mean energy concentration calculation. Figures1, 2, and 3 visualises the comparisons for the interpolated data against the spectrogram distribution set, Choi-Williams distribution set, and Wigner-Ville distribution effect respectively.(a) Turkish earthquake (b) Indonesian earthquake (c) Filipino earthquake (d) New Zealander earthquake

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1: Comparison between interpolated signal of the dataset against heatmap and mean energy concentration of spectrogram/STFT distribution

Table 1 :

 1 Summary of dataset to be processed While the average sampling period can be calculated from observing the summary of each dataset (refer Table

	Earthquake	Range of timestamp	Number of samples Average sampling period
		(UTC)		(seconds)
	Turkish	2023-01-07 09:39:24 to	9399	275.681
		2023-02-06 09:24:52
	Indonesian	2018-09-01 00:25:13 to	2709	924.244
		2018-09-29 23:54:50
	Filipino	2022-07-01 00:37:00 to	1801	1485.441
		2022-07-31 23:45:00
	New	2016-11-01 02:00:55 to	7767	332.316
	Zealander	2016-11-30 22:59:13

Table 3

 3

	).

Table 2 :

 2 Summary of mean energy concentration peaks and corresponding earthquakes of the

	Turkish Earthquake.

Table 3 :

 3 Summary of mean energy concentration peaks and corresponding earthquakes of the

	Indonesian Earthquake.

Table 4 :

 4 Summary of mean energy concentration peaks and corresponding earthquakes of the

	Filipino Earthquake.

Table 5 :

 5 Summary of mean energy concentration peaks and corresponding earthquakes of the New Zealander Earthquake.

Table 6 :

 6 Summary of raw earthquake signal collected from datasets.

	Percentage in
	Turkish Earthquake (%)	5.47 26.1 42.22 13.53 3.19 8.52 0.85 0.09 0.02
	Percentage in
	Indonesian Earthquake (%)	0	0	0.6 19.32 49.7 27.36 2.92	0	0.1
	Percentage in
	Filipino Earthquake (%)	0	0 27.23 53.49 16.36 2.84 0.06 0.06	0
	Percentage in
	New Zealander Earthquake (%)	0 0.45 10.51 61.58 21.74 4.92 0.72 0.06 0.01


Ludovic Goudenège
email: ludovic.goudenege@math.cnrs.fr

Andrea Molent
email: andrea.molent@uniud.it

Antonino Zanette
email: antonino.zanette@uniud.it

Backward Hedging for American Options with Transaction Costs

Keywords: Hedging, Greeks, Transaction costs, American options, Heston model

In this article, we introduce an algorithm called Backward Hedging, designed for hedging European and American options while considering transaction costs. The optimal strategy is determined by minimizing an appropriate loss function, which is based on either a risk measure or the mean squared error of the hedging strategy at maturity. The proposed algorithm moves backward in time, determining, for each time-step and different market states, the optimal hedging strategy that minimizes the loss function at the time the option is exercised, by assuming that the strategy used in the future for hedging the liability is the one determined at the previous steps of the algorithm. The approach avoids machine learning and instead relies on classic optimization techniques, Monte Carlo simulations, and interpolations on a grid. Comparisons with the Deep Hedging algorithm in various numerical experiments showcase the efficiency and accuracy of the proposed method.

Introduction

The groundbreaking work of Black and Scholes [START_REF] Black | The pricing of options and corporate liabilities[END_REF], followed by Merton [START_REF] Merton | Theory of rational option pricing[END_REF], enabled the objective determination of the price of a derivative product as the monetary amount required by an agent for the exact replication of its payoff. This definition assumes that exact replication is achievable, aiming to eliminate any potential risk from holding a short position on the derivative. To achieve this, the authors make several assumptions, the most apparent being the possibility of continuously rebalancing the hedging portfolio. This assumption, however, is not feasible in practice due to physical constraints, making it necessary to rebalance the portfolio at a set of discrete times. These considerations also apply to advanced models, which incorporate additional risk factors, such as Heston's stochastic volatility model [START_REF] Heston | A closed-form solution for options with stochastic volatility with applications to bond and currency options[END_REF]. The infeasibility of continuous-time hedging causes these models to transform into approximate versions of themselves, which are associated with incomplete markets: exact replication of derivative products is generally not possible.

In real-world applications, practitioners typically use continuous-time models to calculate the hedging strategy, despite the impossibility of continuous-time hedging. This is done to exploit the computational power of such models and, specifically, to derive hedging strategies. Furthermore, they include a suitable margin in the price of the derivative to account for the replication error due to discrete-time hedging based on a continuous-time model.

Another limitation to the use of standard theoretical models in a practical context is the presence of frictions, primarily the transaction costs that arise during the construction and rebalancing of the hedging portfolio. Naturally, these costs, while making the replication procedure more expensive on one hand, also discourage the rebalancing prescribed by the theoretical model on the other. Consequently, it is essential to strike a balance in order to minimize the overall cost of the strategy without significantly compromising its effectiveness.

The first studies related to a pricing model in a imperfect market were developed by Duffie and Richardson [START_REF] Duffie | Mean-variance hedging in continuous time[END_REF]. They introduce the mean-variance hedging problem in continuous time, when the hedger can only trade futures contracts on an asset which is positively correlated with the underlying and so they consider the problem of determining the optimal hedge. Schweizer [START_REF] Schweizer | Mean-variance hedging for general claims[END_REF][START_REF] Schweizer | Variance-optimal hedging in discrete time[END_REF] extends the results of Duffie and Richardson and obtains the first ground-breaking result by solving the hedging problem for a general contingent claim when a mean-squared error loss function is considered. Subsequently, Follmer [START_REF] Föllmer | Quantile hedging[END_REF][START_REF] Föllmer | Efficient hedging: cost versus shortfall risk[END_REF] studies theoretical properties of hedging strategies when the loss to minimize is given by the VaR or the CVaR at maturity, computed from the perspective of the hedger. Other authors focus on the development of numerical algorithms for computing optimal hedging strategies. For example, Potter et al. [START_REF] Potters | Hedged Monte-Carlo: low variance derivative pricing with objective probabilities[END_REF] develop an algorithm based on Monte Carlo simulations to determine a hedging strategy that minimizes the quadratic error of hedging between two rebalancing steps. Pochart and Bouchaud [START_REF] Pochart | Option pricing and hedging with minimum local expected shortfall[END_REF] continue this analysis, by considering minimization of a local risk function such as Var and CVaR, also including transaction costs through a first order approximation. The use of convex risk measures for option pricing has been studied also by other authors. For example, Xu [START_REF] Xu | Risk measure pricing and hedging in incomplete markets[END_REF] investigates from a theoretical perspective the evaluation of financial option by risk measure pricing rules. Also Klöppel and Schweizer [START_REF] Klöppel | Dynamic indifference valuation via convex risk measures[END_REF] develop a theoretical investigation and find interesting results on time-consistency about families of conditional convex risk measures. To the same scope, some authors study asymptotic behaviour of prices and strategies. With this regard, Whalley and Wilmott [START_REF] Whalley | An asymptotic analysis of an optimal hedging model for option pricing with transaction costs[END_REF] develop a model for pricing European options in the presence of rehedging transaction costs by means of asymptotic analysis. Kallsen and Muhle-Karbe [START_REF] Kallsen | Option pricing and hedging with small transaction costs[END_REF] obtain asymptotic formulas for utility indifference prices and hedging strategies in the presence of small transaction costs. Similarly, Bouchard et al. [START_REF] Bouchard | Hedging under an expected loss constraint with small transaction costs[END_REF] obtain explicit formulae for option hedging in a market with proportional transaction costs.

Recently, some authors have proposed to use modern reinforcement machine learning methods to compute an optimal hedging strategy through neural networks. One of the main results in this sense is due to Buehler et al. [START_REF] Buehler | Deep Hedging[END_REF]. They address the problem of hedging under conditions of incomplete markets with transaction costs by proposing an algorithm called Deep Hedging: the optimal strategy is determined by training a complex neural network so that it approximates an optimal hedging strategy in terms of an appropriate risk measure. Their approach is greek-free in that the hedging strategy is not determined from the greeks, but rather as the minimizer of a loss function, that is as a solution to a variational problem. The resulting algorithm for finding the optimal hedging strategy has many merits: it is model free, it can manage transaction costs, it can also include, in the feature set, additional information such as trading signals, news analytics or past hedging decisions. This new algorithm has enjoyed considerable success in the financial landscape and beyond. For example, Horvat et al [START_REF] Horvath | Deep Hedging under rough volatility[END_REF] investigate the hedging of short positions on options under rough volatility models again by using the Deep Hedging algorithm. Furthermore, the Deep Hedging algorithm has been recently employed in the insurance field by Carbonneau [START_REF] Carbonneau | Deep hedging of long-term financial derivatives[END_REF], with the scope of hedging long term insurance products such as variable annuities. Despite its merits, there are also some elements that limit the use of the Deep Hedging algorithm: in order to obtain reasonable results, the calculation time required is generally high as the proposed neural network optimises the whole hedging strategy over the entire historical series passed as input. Training such a complex neural network requires considerable numerical effort, and sometimes many iterations are required to converge to the optimal strategy. Moreover, neural networks behave as a black box, which may meet with resistance from regulators (see e.g. Lin [START_REF] Lin | Artificial intelligence, finance, and the law[END_REF]). The Deep Hedging algorithm, as it stands, is specifically designed for European options and it is not suitable for handling American options. In the context of American options, we would like to mention the work by Becker et al. [START_REF] Becker | Pricing and hedging American-style options with deep learning[END_REF], which employs neural networks for pricing and hedging these options.

In this paper, motivated by the algorithm introduced by Buehler et al. [START_REF] Buehler | Deep Hedging[END_REF], we propose an alternative approach that avoids the use of machine learning and relies on traditional optimization and interpolation techniques. We construct an appropriate discretization of the possible market states by using a grid of values, and our algorithm moves backward in time to compute the optimal strategy from each individual state. Consequently, we term this the Backward Hedging (BH) algorithm. The optimal strategy is determined by calculating the optimal quantities of hedging instruments to include in the portfolio, minimizing the desired loss function when applied to forward Monte Carlo simulations generated from the state under consideration. Our proposed method proves to be efficient and reliable, with easily interpretable results. We stress out that, because of lack of time-consistency for most of risk measures (see e.g. Boda and Filar [START_REF] Boda | Time consistent dynamic risk measures[END_REF], Cheridito and Stadje [START_REF] Cheridito | Time-inconsistency of var and time-consistent alternatives[END_REF], Cui et al. [START_REF] Cui | Time-consistent and self-coordination strategies for multiperiod mean-conditional value-at-risk portfolio selection[END_REF]), the problem solved by the Backward Hedging algorithm is not equivalent to the global one, which is solved by the Deep Hedging algorithm, but it can be considered as an alternative problem that best reflects the behavior followed in practice by hedgers, as well as an approximation of the global problem. Furthermore, the use of a backward approach also allows for the valuation of American options. With respect to American options, it is important to note that the approach proposed by Becker et al. [START_REF] Becker | Pricing and hedging American-style options with deep learning[END_REF] differs from ours, as they rely on the exercise strategy derived from the continuous-time model, rather than adopting the worst-case strategy for the hedger, as we suggest in our proposal.

Our investigation includes numerical tests in the one-dimensional Heston model, primarily considered in [START_REF] Buehler | Deep Hedging[END_REF]. The results obtained are consistent with those deduced from the Deep Hedging algorithm and are competitive in terms of quality, confirming the validity of our proposed procedure. Specifically, although the solution of the backward problem is usually not optimal for the global problem, numerical results show that the hedging strategy returned by our algorithm is effective and sometimes even a better minimizer for the global problem than the one obtained with the Deep Hedging algorithm. It is worth noting that Buehler et al. [START_REF] Buehler | Deep Hedging[END_REF] demonstrate the applicability of the Deep Hedging algorithm in the context of a multidimensional Heston model. Our proposed algorithm, however, cannot efficiently address this case due to the curse of dimensionality, which can only be solved by machine learning-based algorithms. Nevertheless, we believe that our approach may be useful for practitioners when considering low-dimensional problems.

The paper is organized as follows. Section 2 introduces the pricing framework of the Deep Hedging algorithm. Section 3 discuss the proposed new algorithm. Section 4 is devoted to the experimental results. Finally, the conclusions follow in Section 5.

Deep Hedging framework

The hedging model presented in this manuscript is developed within the framework discussed in the Buehler et al. paper [START_REF] Buehler | Deep Hedging[END_REF]. The Deep Hedging algorithm is also used as a point of comparison for our algorithm; therefore, we begin our discussion by briefly presenting it. Although the algorithm is highly adaptable to different stochastic models, we present it with reference to Heston's (one-dimensional) stochastic model, which is also studied in the original manuscript. It is worth noting that in the original version of the algorithm, the authors assume both the drift rate and the interest rate to be zero. Here, we present an updated version of the Deep Hedging algorithm, which works for non-zero drift and interest rate.

The Deep Hedging algorithm allows one to approximate optimal hedging strategies for managing a short position on a European option. The search for an optimal strategy is carried out by exploiting neural networks and it is based on the minimization of a certain loss function. The optimal strategy ultimately determines the option price as the amount of money required to make such a strategy acceptable according to a risk profile. This model makes it possible to determine the value of an option even when perfect hedging is not possible, also when trading costs are considered.

Settings

We consider a financial market with finite time horizon T and trading dates 0 = t 0 < t 1 < . . . < t N = T , with t n = n • ∆t for n = 0, . . . , N . Let S 1 be the stochastic process representing the price of a stock and V the associated variance process. These two stochastic processes are supposed to evolve, under a real-world probability P, according to the Heston model, that is the following stochastic differential equations hold:

S 1 t = µS 1 t dt + √ V t S 1 t dB t , S 1 0 = s 0 , V t = α (b -V t ) dt + σ √ V t dW t , V 0 = v 0 ,
where B and W are one-dimensional Brownian motions with their correlation ρ BW a real constant, just as µ, α, b, σ, s 0 and v 0 . In the original version of the manuscript [START_REF] Buehler | Deep Hedging[END_REF], the authors consider a market with a zero drift and interest rate, that is µ = r = 0, however, here, following Rémillard and Rubenthaler [START_REF] Rémillard | Optimal hedging in discrete and continuous time[END_REF], we allow these values to be non-null. Unlike process S 1 , process V is not directly tradable, nonetheless we assume that an idealized variance swap is available in the market, and we denote its price by S 2 . In particular, the maturity of the variance swap is the same as the considered derivative, i.e. T , and if we define the filtration

F t := σ S 1 s , V s : s ∈ [0, t] , then S 2 t := e -r(T -t) • E Q T 0 V s ds | F t , t ∈ [0, T ].
with Q a risk-neutral probability. Moreover, if we assume the market price of risk for the variance process to be zero, one can prove that the idealized price of the variance swap writes down

S 2 t = e -r(T -t) • t 0 V s ds + L (t, V t )
where the function L is defined by

L(t, v) = v -b α 1 -e -α(T -t) + b(T -t).
Then, S 1 and S 2 are the prices of tradeable assets that can be used to set up a portfolio to hedge a short position on a derivative. Now, let us consider a European option, whose payoff is a random variable F T measurable, denoted by φ S 1 T . We assume the position of an economic agent who sold the option, so φ S 1 T is a liability, and we have to hedge such a short position.

Trading strategies

In order to hedge a liability φ S 1

T at T , one trades S 1 and S 2 by using an R 2 -valued F-adapted stochastic process δ = (δ n ) n=0,...,N -1 with δ n = δ 1 n , δ 2 n . Please observe that, δ i n denotes the agent's holdings of the i-th asset, after rebalancing the portfolio at time t n . For simplicity of notation, we may also define δ i -1 = δ i N := 0 for i = 1, 2, which implies an empty portfolio at beginning and full liquidation at maturity.

Let p 0 denote the cost for building the portfolio at t 0 , so p 0 is the amount of money required in payment by the buyer of the option. In a market without transaction costs the agent's global profit and loss P L T 0 at time T is thus given by

P L T 0 φ S 1 T , p 0 , δ = -φ S 1 T + p 0 e rT + (δ • S) N 0 ,
where

(δ • S) n2 n1 := n2-1 k=n1 δ 1 k S 1 t k+1 -S 1 t k e r∆t + δ 2 k S 2 t k+1 -S 2 t k e r∆t e (T -t k+1 )r
is the profit due to the hedging strategy, with respect to the time intervals [t k , t k+1 ] for k = n 1 , . . . , n 2 -1, capitalised to maturity.

Trading costs can also be included in the model. We assume that any trading activity causes proportional costs as follows: if the hedger decides to buy or sell α ∈ R units of S i at time t k , then he incurs in a transaction cost which amounts to ε • S i k • |α|, with ε the proportional cost for trading 1 currency unit of asset. The capitalised to maturity cost of the trading strategy δ, with respect to trading times from t n1 to t n2 , is therefore given by

C n2 n1 (δ) := n2 k=n1 ε S 1 t k δ 1 k -δ 1 k-1 + S 2 t k δ 2 k -δ 2 k-1 e (T -t k )r .
We stress out that, according to the original paper, C N 0 (δ) includes the cost of setting up the portfolio at time 0 and the costs of divestment at time T . As a result, when transaction costs are applied, the agent's terminal portfolio value at T reads out

P L T 0 φ S 1 T , p 0 , δ := -φ S 1 T + p 0 e rT + (δ • S) N 0 -C N 0 (δ).
Such a random variable can be considered as the error of the trading strategy, and one seeks to minimize it.

Hedging and pricing

The basic situation to considered is given by a complete market with continuous-time trading, no transaction costs, and unconstrained hedging. In this particular case, for any liability φ S 1 T , there exists a unique replication strategy δ and a unique fair price p 0 ∈ R such that

-φ S 1
T + p 0 e rT + (δ • S) N 0 = 0 holds P-a.s.. More generally, the optimal hedging strategy δ and the option price p 0 can be defined as the solution of the following problem: computing p 0 ∈ R and δ ∈ H such that

E -φ S 1 T + p 0 e rT + (δ • S) N 0 2 = 0, (2.1)
where H stands for the set of all possible trading strategies.

In an incomplete market with frictions, perfect hedging is (usually) no more possible, so, it is no longer possible to obtain a strategy that generates an almost certainly zero error. Consequently, the new objective is to minimize the error of the hedging strategy, so, one considers the problem of minimizing the mean squared error of the strategy:

inf p0∈R inf δ∈H E -φ S 1 T + p 0 e rT + (δ • S) N 0 -C N 0 (δ) 2 . (2.2)
The solution p 0 of (2.2) can be regarded as the fair hedging price, but it does not take into account the non-eliminable risk faced by the hedger. In this case, one has to specify a certain criterion that defines an acceptable "minimal price" for any liability, a price that includes compensation for the risk assumed by the hedger. Such a price, can be defined as the minimum amount of money to add to the hedging portfolio at the launch of the contract, so that, when the optimal strategy is employed, the overall position becomes acceptable at maturity. The criterion for identifying what is an acceptable position is defined by means of a convex risk measure ρ, such that ρ (0) = 0. For example, hereafter, we consider ρ as the α-level CVaR, that is ρ (X) = CVaR α (X). The problem to solve is then computing

inf δ∈H ρ -φ S 1 T + (δ • S) N 0 -C N 0 (δ) , (2.3)
so that the hedger's price writes down

p0 = e -rT inf δ∈H ρ -φ S 1 T + (δ • S) N 0 -C N 0 (δ) . (2.4)
Let δ ∈ H be the optimal strategy, solution of problem (2.2) or (2.3), according to the evaluation approach under investigation. Please, observe that the optimal hedging strategy δi n at time t n can be seen as the outputs of a function that depends on filtration at time t n : δi n = δi n S 1 0 , S 2 0 , . . . , S 1 n , S 2 n . Accordingly, problems (2.2) and (2.3) can be interpreted as variational problems, which can be tackled by exploiting a reinforcement learning algorithm based on neural networks.

The Deep Hedging algorithm

Neural networks can be used to solve variational problems such as (2.2) and (2.3) by choosing a suitable loss function.

The first step consists in simulating the paths of the underlyings S 1 , S 2 at discrete times t n . These simulations are used to train a semi-recurrent neural network: at each time-step t n , for n = 0, . . . , N -1 it takes as inputs the 4 values of S 1 tn , S 2 tn , δ 1 n-1 and δ 2 n-1 and it outputs δ 1 n and δ 2 n . The network is trained to minimize the loss function, which value is determined based on the simulations of the paths. At the end of this process, the network returns a hedging strategy, denoted by δ, that approximates the solution to the problems (2.2) or (2.3). Specifically, the weight matrices and biases of the neurons are the parameters of the neural network but also the initial values δ 1 0 and δ 2 0 are added among the parameters to be optimize. Now, let θ be the vector all the parameters of the neural network. So, θ is computed as a minimum point of a suitable loss function ℓ (θ). While problem (2.3) is considered, the loss function is

ℓ (θ) = ρ -φ S 1 T + (δ • S) N 0 -C N 0 (δ) .
As far as problem (2.2) is considered, also p 0 is included among the parameters to be optimized and the loss function is set as

ℓ (θ) = E -φ S 1 T + p 0 e rT + (δ • S) N 0 -C N 0 2 .
The neural network is then trained by using the ADAM algorithm (Kingma and Ba [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF]). Finally, when considering problem (2.3), the trained neural network is tested on out-of-sample dataset to determine the price p 0 , so that the measurement of network performance is not altered by overfitting.

Once the neural network has been trained, it returns the following outputs: the option prince p 0 and the initial values δ1 0 and δ2 0 as well as it can be employed to calculate delta values at each possible trading timestep. This means that the neural network can be used for hedging purposes: given a path S 1 tn , S 2 tn n=0,...,N or even a partial path S 1 t k , S 2 t k k=0,...,n (with n ≤ N ), the neural network returns its optimal values for δ1 k , δ2 k k=0,...,n , so one can effectively exploit the Deep Hedging algorithm for hedging a short position.

Backward Hedging

The Deep Hedging algorithm relies on neural networks to achieve global optimisation of the hedging strategy. This means that the algorithm simultaneously optimises all the parameters θ that define the entire hedging strategy. This approach results in a single high-dimensional optimisation problem.

The algorithm we propose, on the other hand, solves a large number of optimisation problems in only two variables, that is δ 1 n , δ 2 n . Starting from a specific configuration of the market state, it determines the optimal delta values for that specific configuration. Thanks to linear interpolation, the hedging strategy is extended to the entire space of interest. We emphasise that, unlike other works already present in the literature (see e.g. Potter et al. [START_REF] Potters | Hedged Monte-Carlo: low variance derivative pricing with objective probabilities[END_REF]) that consider the loss between two trading times, the loss function to be optimised is evaluated at maturity, i.e. it considers the entire hedging strategy, from the considered time-step until maturity.

Finally, we would like to point out that the pseudocodes of the main versions of the Backward Hedging algorithm, discussed below, are included in the concluding Appendix A.

The backward approach

We begin our discussion by presenting the proposed hedging model for the special case ε = 0, that is no transaction costs and loss function based on the risk measure ρ, as in problem (2.3).

The hedging strategy we propose is determined by proceeding backward in time, starting from time t N -1 , to the initial instant t=0. At the generic hedging instant t n , the amount of hedging instruments to be held in the portfolio δ 1 n and δ 2 n are determined by minimizing an appropriate measure of risk relative to the future hedging error, that is, between t n and T . Specifically, the Backward Hedging strategy δ at time t n is defined as the solution of the following problem:

δ1 n , δ2 n = arg min (x,y)∈R 2 ρ -φ S 1 T + ( δ • S) N n+1 + x S 1 tn+1 -S 1 tn e r∆t + y S 2 tn+1 -S 2 tn e r∆t e (T -tn+1)r . (3.1)
Let us remark that in equation (3.1), at time t n , S 1 tn and S 2 tn are measurable, while S 1 t k and S 2 t k for k > n are random variables. We, observe that problem (3.1) does not involve the previous hedging strategy, but only the future hedging strategy, i.e. δ1 k , δ2 k for k = n+1, . . . , N . We stress out that excluding past history seems reasonable from a practical perspective since an hedger decides his strategy looking only at future results. In fact, the gain and the losses already realised through the hedging strategy are measurable values and, due to the cash invariance property of risk measures, their inclusion/exclusion does not change the solution to problem (3.1). Moreover, problem (3.1) is suitable to be tackled through a backward procedure: if one solves this kind of problems backward in time, future hedging strategy should be available as referred to time steps already treated. Furthermore, the applicability of a backward approach is also ensured by the fact that Heston's model is Markovian, so the future evolution of processes, conditioned by the state of the present market, is independent of the past. We can therefore observe that values δ1 k , δ2 k only depend on the current state of the market. Once one has solved problem (3.1) for t n = t N , . . . , 0, the Backward Hedging option price can be computed as

p0 = e -rT ρ -φ S 1 T + ( δ • S) N 0 -C N 0 ( δ) . (3.2)
Before presenting the resolution of the problem (3.1) in Heston's model, we discuss the relation between the strategies δ and δ, that is the solution of the global hedging problem (2.3) and of the backward problem (3.1), respectively. We now show that the two strategies, in general, are not equivalent.

A key example

We consider, as an example, a discrete time model, with only 3 trading times, namely t 0 = 0, t 1 = 1 and t 2 = 2. For this example only, the market includes a single underlying with initial value S 0 = 1. The underlying can assume 5 equiprobable values at time t 1 , which are S i 1 = S 0 • u 2i-6 , for i = 1, . . . , 5 and u = exp 0.2 1/4 . Moreover, for each value of S i 1 , there are 5 equiprobable continuation values, namely S i,j 2 = S i 1 • u 2j-6 , for j = 1, . . . , 5. So there are 5 possible paths from t 0 to t 1 and 25 possible paths from t 0 to t 2 . We stress out that these values are the same as those generated at time-steps 0, 4 and 8 by the CRR model with T = 2 and N = 8 discretization steps, with the only difference being that, to simplify the calculations, we assume transaction probabilities between two consecutive instants all equal to 1/5. Also, for simplicity, we assume the interest rate r to be zero. In this model, an hedging strategy consists of 6 values: δ 0 , i.e. the amount of stocks to include in the hedging portfolio at time t = 0, and δ 1,i = δ 1 (S i 1 ), i = 1, . . . , 5, i.e. the amount of stocks to include in the hedging portfolio at time t 1 according to the 5 possible values of the underlying at that time. Now, we consider a short position on a strangle, i.e. a short position on a call and on a put, with strikes K C = 0.8 and K P = 1.3 respectively. In our example, we want to hedge the strangle so as to minimize the 40%-CVaR. For the global-optimal strategy δ, this consists of choosing the 6 hedging parameters so as to minimize the CVaR relative to profit-and-loss between t 0 and t 2 , which we denote by CVaR 2 0 . For the backward strategy δ, on the other hand, this consists of determining at time t 1 , for each of the possible values of the underlying S i 1 , the hedging parameter δ1 (S i 1 ) so as to minimize the CVaR relative to the profitand-loss between t 1 and t 2 , which we denote by CVaR 2 1 , and then the hedging parameter δ0 at time t 0 so as to minimize the CVaR 2 0 between t 0 and t 2 , being subject to the constraint of using δ1 (S strategy at time t 1 . Let us point that, since we are assuming equal transaction probabilities and 40% as CVaR level, then CVaR 2 0 is equal to the opposite of the average of the 10 lower profit-and-loss outcomes between t 0 and t 2 , while CVaR 2 1 is equal to the opposite of the average of the 2 lower profit-and-loss outcomes between t 1 and t 2 . The two optimal strategies are shown in Table 1, while Table 2 reports the computations for those strategies. Specifically, the first three columns show the 5 possible values of S 1 , the 25 final values of S 2 and their payoffs, respectively. The next columns report the profit and loss values for each of the 25 possible trajectories, calculated by applying the 2 strategies under consideration. These profit-and-loss are calculated between t 0 and t 2 , denoted by P L 2 0 , or between t 1 and t 2 , denoted by P L 2 1 . Values in bold are the worst values, and they contribute to the calculation of CVaR values. For both the two strategies, CVaR 2 1 is reported for each of the 5 possible values of S 1 , while the two values for CVaR 2 0 are reported at the bottom of the Table. We can observe that the two strategies are successful in their respective problems: δ generates a smaller CVaR 2 0 than δ, while δ generates smaller (or equal) CVaR 2 1 for all possible values of S 1 . We can therefore conclude that the two strategies δ and δ, and therefore the problems (2.3) and (3.1), are not equivalent.

Comparison between strategies

At this point, it is worth asking which of the two strategies is the most correct to consider for practical purposes. The global optimal strategy, pursued by the Deep-Hedging algorithm, hides a pitfall. The strategy δ is optimal at time zero, that is between t 0 = 0 and T , but only at time zero and has no optimality property at times t n > 0, between t n > 0 and T : in fact, as shown by the previous example, the strategy δ is not optimal at time t 1 , as it does not always minimize CVaR 2 1 , in as exceeded by δ. If we assume that, at any time t n > 0, the hedger wishes to minimize the CVaR between t n and T , then he will be led to abandon strategy δ that he had decided at time 0, in favor of a more performing one. The initial valuation is therefore not consistent with what is expected of the hedger: the option is priced based on a strategy which will not be followed afterwards. This annoying phenomenon is related to the absence of time-consistency for CVaR, and more generally for the most common risk measures. The δ strategy, on the other hand, at time t n is optimal given the strategy that will be followed subsequently, therefore it is not convenient to modify it, unlike δ, and it is less penalizing in future times. Once established, it should be maintained (unless it is replaced with another δ strategy recomputed to be optimal between t n and T , which however will be disregarded in subsequent times t k > t n ). This aspect could be appreciated by the supervisory bodies which must validate the hedging procedure. Finally, we point out that Basak and Chabakauri [START_REF] Basak | Dynamic mean-variance asset allocation[END_REF] and subsequently Cui et al. [START_REF] Cui | Time-consistent and self-coordination strategies for multiperiod mean-conditional value-at-risk portfolio selection[END_REF] considerations in the context of portfolio allocation: although the global-optimal strategy first appears to have strong global investment performance, investors would deviate from such a strategy throughout the investment process.

For the reasons previously mentioned, we propose to use δ in place of δ .

European option with no transaction costs, based on risk measure

We discuss the numerical procedure we propose to solve problem (3.1). The issues that one has to face are controlling the hedging strategy for future time-steps, computing the objective function, that is the risk measure, and optimizing the objective function.

First of all, we propose to estimate the risk measure ρ through (forward) M Monte Carlo simulations, from t n to T , of the underlying processes S 1 and S 2 conditioned by F tn . In order to simulate future values of S 1 , S 2 , one only needs to know S 1 tn , V tn and I tn = tn 0 V s ds in place of all filtration F tn , so the market state can be described by this triple. Please, note that knowing S 1 tn , S 2 tn is not sufficient to simulate future values S 1 t k , S 2 t k |F tn , k = n + 1, . . . , N , as, to simulate the underlying S 1 t k , one must first simulate the volatility V t k , and the initial value of the volatility V tn cannot be deduced from S 1 tn , S 2 tn . Therefore we also generate the values of I and V , and denote the simulated values with

S n = S 1 k,m , S 2 k,m , I k,m , V k,m , k = n + 1, . . . , N, m = 1, . . . , M .
Now, for each time-step t n , we consider a grid of values

G n = G S n × G I n × G V n which is defined as follows: G S n
is a uniform mesh of N S points for S 1 tn , G I n is a uniform mesh of N I points for I tn and G V n is a mesh of N V points for V tn . The maximum and minimum values of G S n and G I n are set by computing appropriate quantiles (in our tests we consider 0.1% and 99.9% quantiles) for S 1 tn and I tn . Those quantiles are estimated from M 0 Monte Carlo simulations of the Heston model, starting from the initial values S 1 0 = s 0 and V 0 = v 0 , which we denote with Ŝn = Ŝ1 k,m , Îk,m , k = 1, . . . , N, m = 1, . . . , M 0 . Simulations of the stochastic processes are obtained through the hybrid Monte Carlo-Tree method introduced by Briani et al. [START_REF] Briani | A hybrid approach for the implementation of the Heston model[END_REF][START_REF] Briani | A hybrid tree/finite-difference approach for Heston-Hull-White-type models[END_REF], together with the trinomial tree introduced by Goudenège et al. [START_REF] Goudenège | Gaussian process regression for pricing variable annuities with stochastic volatility and interest rate[END_REF], by using N T simulation time-steps (N T must be a multiple of N , so that the simulated values at the trading epochs can be obtained by subsampling). This simulation method is particularly interesting in that the possible states for the random variable V tn are discretized through the use of a trinomial tree T that exactly reproduces the first and second moments of the transition distribution of the process V , thus leading to accurate results with only a few time-discretization steps. This trinomial tree is also used to define the grid G V n as the values of V associated to the tree nodes at time t n . Moreover, it is worth observing that the simulated values of V always belong to the grid. Finally, we also set G 0 = {s 0 } × {0} × {v 0 } for the market state at inception. Problem (3.1) is solved for each point of G backward in time so, at time t n , the hedging strategy is known at future time-steps {t k } k=n+1,...,N -1 , but only for those values of S 1 t k , I t k , V t k that are points of G k . However, the simulated values of S t k and I t k for the estimation of the risk measure may not belong to the grid G k and therefore the hedging strategy would not be known for those special values. This problem is overcome by using linear interpolation, which allows one to extend the hedging strategy from G k to all R 3 . In this regard, we point out that it is necessary to resort to an interpolation in the S and I coordinates only, as the simulated values of V belong to the grid thanks to the use of the trinomial tree T for simulation. We also point out that although other forms of interpolation are possible (see, e.g., Akima [START_REF] Akima | A method of bivariate interpolation and smooth surface fitting based on local procedures[END_REF][START_REF] Akima | A method of bivariate interpolation and smooth surface fitting for irregularly distributed data points[END_REF] for Modified Akima interpolation or Press et al. [START_REF] Press | Numerical recipes 3rd edition: The art of scientific computing[END_REF] for Bicubic Interpolation), linear interpolation is chosen for its low computational cost. This element is decisive as the proposed algorithm requires the calculation of a large number of interpolations. As a result, we can represent the hedging strategy as a collection of function δ1 n , δ2 n : R 3 → R 2 such that δ1 n , δ2 n (s, i, v) are the optimal hedging values when S 1 tn = s, I tn = i and V tn = v. Now that we have clarified how to calculate the objective function, we must solve the optimization problem (3.1), thus determining the specific values δ1 n , δ2 n for the grid point under evaluation. This is done by exploiting a standard numerical procedure. Specifically, we used the Simplex Search method of Lagarias et al. [START_REF] Lagarias | Convergence properties of the Nelder-Mead simplex method in low dimensions[END_REF], which is an efficient derivative-free algorithm and it is particularly suitable for problems of this kind.

Summing up, the Backward Hedging algorithm proceeds backward in time by solving problems (3.1) for all nodes in the grid, starting from time t N -1 down to time t 0 . Once the hedging strategy has been determined for all trading times, we simulate a new sample of M 1 Monte Carlo simulations

S0 = S1 k,m , S2 k,m , Īk,m , Vk,m , k = 1, . . . , N, m = 1, . . . , M 1
and we apply the algorithm on this new out-of-sample data. The risk measure determined on this new sample defines the initial option price p0 with respect to the Backward Hedging algorithm.

Loss as mean squared error

The Backward Hedging algorithm can also be employed to tackle problem (2.2), by a backward approach. Again, let us again assume that no transaction costs are applied. It is easy to prove that the optimal value of p 0 according to (2.2) is given by

p 0 = -e -rT E -φ S 1 T + (δ • S) N 0 , (3.3)
and for this choice of p 0 , we obtain

E -φ S 1 T + p 0 e rT + (δ • S) N 0 2 = Var -φ S 1 T + (δ • S) N 0 , thus (2.2) reduces to solving inf δ∈H Var -φ S 1 T + (δ • S) N 0 . (3.4)
Now, let us consider the backward version of problem (3.4), that is at time t n one has to solve the following problem

δ1 n , δ2 n = arg min (x,y)∈R 2 Var -φ S 1 T + ( δ • S) N n+1 + x S 1 tn+1 -S 1 tn e r∆t + y S 2 tn+1 -S 2
tn e r∆t e (T -tn+1)r .

(3.5) Problem (3.5) is of the same type of (3.1), with the only difference being that ρ is replaced with Var, so one can tackle it again with the Backward Hedging algorithm, as described in the previous Subsection.

American options

American options can be exercised at any moment before maturity. With the scope of pricing, we approximate such a feature with a weaker one: we suppose that the option can be exercised at any trading time. If the option is exercised, the hedger pays the payoff to the buyer, otherwise he rebalances the hedging portfolio. Pricing requires to model the exercising strategy, which depends on the buyer. Here, we assume the so called worst-hedging-case, which means that the option is exercised in the worst case for the hedger, that is as soon as the payoff is higher than the hedging cost, which corresponds to the continuation cost. Specifically, the option price pn at time t n is define as a function of S 1 tn , I tn , V tn , again as the cost for hedging the position, from t n up the the exercise time τ n that is

pn S 1 tn , I tn , V tn = inf δ∈H ρ -φ S 1 τn + ( δ • S)
τn/∆t 0 e -r(τn-tn) .

(3.6) Specifically, τ n is the stopping time defined as

τ n = min τ ∈ {t n , . . . , T } s.t. φ S 1 τ ≥ pn S 1 τ , I τ , V τ . (3.7)
and so we can tackle such a problem again by backward induction. In particular, the hedging strategy is computed by solving the following problem

δ1 n , δ2 n = arg min (x,y)∈R 2 ρ -φ S 1 τn + ( δ • S) τn/∆t n+1 +
x S 1 tn+1 -S 1 tn e r∆t + y S 2 tn+1 -S 2 tn e r∆t e r(τn-tn+1) . (3.8)

Again, such a problem is solved backward in time, by exploiting Monte Carlo simulations and linear interpolation. In particular, the calculation of the exercise strategy requires the comparison between the payoff and option price, so, in this case, linear interpolation is used for both interpolating the hedging strategy and the option price at each trading time.

Transaction costs

Let us consider proportional transaction on traded assets as in Buehler et al. [START_REF] Buehler | Deep Hedging[END_REF]. The inclusion of transaction costs complicates the evaluation of the hedging strategy because the optimal strategy at one instant is closely related to the strategy employed at the earlier time-step. We start by considering the problem of hedging a European derivative, based on the risk measure ρ, when proportional transaction costs are due. The following equation defines the Backward Hedging strategy at time t n :

δ1 n , δ2 n = arg min (x,y)∈R 2 ρ ( -φ S 1 T + ( δ • S) N n+1 + C N n+1 ( δ) (3.9)
+ x S 1 tn+1 -S 1 tn e r∆t + y S 2 tn+1 -S 2 tn e r∆t e (T -tn+1)r +ε S 1 tn x -δ 1 n-1 + S 2 tn y -δ 2 n-1 e (T -tn)r .
Similarly to what we did for the no-transaction-costs case, we replace conditioning with respect to F tn with a discrete set of data, i.e. S 1 tn = s,

I tn = i, V tn = v, δ 1 n-1 = d 1 , δ 2 n-1 = d 2 .
We stress out that, in this case, we cannot eliminate the formula's dependence on δ 1 n-1 , δ 2 n-1 , which, proceeding backward, will be computed at the next time-step as δ1 n-1 , δ2 n-1 , with the only exception for δ 1 0 = δ 1 0 = 0. So, we cannot apply the backward approach directly. To overcome this obstacle, we consider a 5 dimensional grid, which discretizes the conditioning event. Specifically, for each time-step t n , we consider a grid of values

G n = G S n × G I n × G V n × G δ 1 n × G δ 2
n which is defined as follows: G S n , G I n and G V n are defined as in Section 3.4, while G δ 1 n and G δ 2 n are two uniform mesh of N δ 1 and N δ 2 points that represent the values for δ 1 n-1 and δ 2 n-1 respectively. Specifically, the maximum and minimum values of G δ 1 n and G δ 2 n are determined as specific quantiles (in our tests we consider 0.1% and 99.9% quantiles) of δ 1 n-1 and δ 2 n-1 when ε = 0. Basically, those quantiles are computed by exploiting Monte Carlo simulations of the Backward Hedging algorithm for the case without transaction costs, starting from time t = 0. The interval for the grids G δ 1 n and G δ 2 n are appropriately extended to include the zero value, which is the delta value at the beginning and at the end of the hedging strategy.

The Backward Hedging algorithm for this case is developed following the same valuation principle adopted for the case without transaction costs: it proceeds backward in time and, at each time-step t n , for each point of the 5-dimensional, it determines the optimal values of δ 1 n and δ 2 n through a series of Monte Carlo simulations up to maturity (or to the exercise time with respect to American options), considering the values associated with the starting grid point as the market state parameters. Again, linear interpolation is used to extend the hedging strategy from the grid nodes to the entire R 5 space. Furthermore, in order to favour the stability of the algorithm in the case of extreme trajectories (that is simulated values are not inside in the grid), the values predicted by the interpolation for δ 1 k and δ 2 k for k = n + 1, . . . , N , are restricted to the discretization interval of G δ 1 k and G δ 2 k respectively. Remark 1. In all cases considered in the preceding sections, for each time-step, the calculation of the hedging strategies for each grid point are independent of each other, so this step can be parallelised, thus reducing the total calculation time of the Backward Hedging algorithm.

Numerical results

In this Section we discuss some numerical results related to the proposed algorithm. Specifically, we compute the optimal hedging strategy for a European or American put option, considering either the minimization of the mean square error of hedging or the CVaR at maturity as the objective, including (i.e. ε = 1%) or not transaction costs (i.e. ε = 0). Moreover, we consider different settings for N , the number of trading dates, and M , the number of Monte Carlo simulations used for solving the optimization problems that define the hedging strategy. In particular, in order to obtain reasonable calculation times, we have considered values M = 10 3 , 10 4 , 10 5 . We first consider the model parameters employed for a test by Buehler et al. [START_REF] Buehler | Deep Hedging[END_REF], and secondly the same parameters but with positive interest rate and drift. The employed values are reported in Table 3 together with the put option parameters, while Table 4 displays the parameters for the Backward Hedging algorithm.

The following tables show the hedging prices at time t 0 in the various cases considered, obtained by means of the Backward Hedging algorithm (BH), which has been implemented in the MATLAB environment. For the sake of comparison, these prices were also calculated using our implementation of the Deep Hedging algorithm (DH), which has been implemented in Python by us. The parameters for the Deep Hedging algorithm are listed in Table 5. We point out that, for the Deep Hedging algorithm, there are two parameters that govern the number of Monte Carlo simulations: M the number of simulations for training the neural network and M 1 the number of out-of-sample simulations performed at the end to generate the option price. For the Backward Hedging algorithm, to these two parameters, we add M 0 , the number of simulations used to determine the grid that discretizes the market states. For the sake of comparison, we also report the computational times, although these times should be taken with some detachment as they refer to codes developed in different environments. Furthermore, the implementation of the Deep Hedging algorithm was carried out by us and there is probably room for improvement.

The timed comparisons were carried out on a personal computer equipped with an Intel i5-1035G1 processor and 8 GB of RAM, and only one core was used in order to stabilize the calculation times as much as possible.

We emphasise that the prices shown in the following tables are determined by considering the profit and loss of the hedging strategy from time zero to maturity, thus representing an opportunity to compare the two hedging strategies on the same objective functions, ρ P L T 0 and Var P L T 0 . In this regard, DH's strategy should be the optimum point in the set of hedging strategies for the global problem, while BH's strategy can be considered as a solution to the same problem, which includes the constraint of given by the recursive scheme that defines the backward strategy. Then, through these numerical tests, we can compare the two strategies and see how optimal they are for the global problem.

Symbol

Mean squared error

In this Section we consider the minimization of the mean squered hedging error as target. Tables 6 and7 show the results for the European and the American put option for r = µ = 0, while Tables 8 and9 report the corresponding results for the r = µ = 0.1 case. We also show the price in Heston's continuous model, obtained through the Hybrid Tree-PDE algorithm by Briani et al. [START_REF] Briani | A hybrid approach for the implementation of the Heston model[END_REF][START_REF] Briani | A hybrid tree/finite-difference approach for Heston-Hull-White-type models[END_REF]. This price serves as a benchmark for the case where transaction costs are zero. No benchmark is available for positive transaction costs. In all the considered cases, the prices returned by th BH algorithm are very stable with respect to M and the results are consistent with both the benchmarks (when available) and the results returned by the DH algorithm (for the European option). When r = µ = 0, the American price is very similar to the European price, but when r = µ = 0.10, the difference is clear. In the latter case, we observe that the price increases as the number of instants of hedging increases because the number of possibilities for exercising the option increases.

To conclude, the proposed algorithm proves to be very efficient and stable.

ε = 0 ε = 1% Benchmark 3.7165 - M N 4

CVaR error

In this Section we consider the minimization of the 95%-CVaR of the hedging strategy at maturity as target. Before discussing the results, we point out that, in order to facilitate comparison between the proposed algorithms, in this case we have also included, for the DH algorithm only, the case M = 10 6 , which is necessary to obtain competitive results with BH. Tables 10 and11 show the results for the European and the American put option for r = µ = 0, while Tables 12 and13 report the corresponding results for the r = µ = 0.1 case. No benchmark is available in this case. In all cases considered, the prices returned by th BH algorithm are very stable with respect to M . As far as European options are considered, the price decreases as N , the number of possible portfolio rebalances, increases, since the higher this number, the smaller the error committed at the end of the hedging period and the lower the relative CVaR. In fact, when one aims to minimise a risk measure of replication error, unlike mean squares error minimisation, the lower the price obtained, the better the strategy. As observed in the previous numerical experiments, discussed in Subsection 4.1, in the presence of transaction costs, the option price increases, as expected. As noted in the previous section, when r = µ = 0, the American price is very similar to the European price, but when r = µ = 0.10, the difference is clear. In the latter case, we observe that the price decreases: while the increase in trading dates favours the optimal exercise of the option that raises the price, it also favours the development of an effective hedging strategy that lowers the price. Between these two contrasting effects, the second seems to prevail. As shown in Tables 10 and12, the comparison with DH shows that BH's proposed strategy is also very efficient in dealing with the global problem: for all the parameter configurations considered, the price of the BH strategy is very similar to the price of DH's strategy, as their difference is small, and in some cases even lower. Moreover, the computational time tends to reward BH. We also observe that the run times of BH are more regular than those of DH, as the stop criterion employed for training the neural network seems to be more sensitive to specific input data. In particular, DH times seem to increase up to M = 10 5 , whereas for M = 10 6 they are of the same order of magnitude as for M = 10 5 . In conclusion, BH again proves to be very efficient and stable.

To complete our analysis, we present Figure 4.1 which shows a plot of the optimal values of δ1 3 and δ2 3 , the hedging strategy at time t 3 = 0.75, computed through the Backward Hedging algorithm, as a function of the hedging strategy at time t 2 = 0.5, δ 1 2 , δ 2 2 . In particular, we consider an American put option with N = 4 hedging steps, S 0.75 = 100, I 0.75 = 0.04 • 0.75 and V 0.75 = 0.04, with M = 10 4 Monte Carlo simulations. Hedging is based on the CVaR minimization and it employs ε = 1% proportional transaction cost. To improve the graphic detail, the values N δ1 and N δ2 have been set equal to 11.

The graphs suggest that both δ1 3 and δ2 3 increase as δ 1 2 and δ 2 2 increase, with a greater dependence on δ 1 2 . This fact probably results from the low absolute cost in trading S 2 due to the small price of S 2 with respect of S 1 . To confirm this, in Figure 4.2 we can observe a graph similar to that of Figure 4.1, in which the relative cost of trading S 2 has been raised to 10%, while that for S 1 has remained 1%. In this case, two increasing graphs are still observed, but the dependence on δ 1 2 is more relevant, in particular for δ2 3 . Table 13: Hedging cost for an American put option with r = µ = 0.1, while minimizing loss based on CVaR. . Proportional cost for trading S 1 is 1%, while the proportional cost for trading S 2 has been raised to 10%.

ε = 0 ε = 1% M N 4

Conclusions

In this paper, we have discussed a method for identifying the most efficient hedging strategies for derivative products, including both European and American options, while taking into account transaction costs. The hedging strategy, defined as the quantities of assets held in the hedging portfolio, is determined to minimize a risk measure or the mean squared error of the hedging error at the time the option is exercised. We calculate the optimal strategy for each point on a grid by using Monte Carlo simulations and then we extend the strategy to other points via interpolation. Although the proposed strategy is not optimal at time zero for the global optimization problem, it offers, at each time-step, a strategy that is optimal given the strategy that will be followed later, so it is easier to follow in practice, pleasing both practitioners and supervisors. The proposed method proves to be effective and reliable, offering clear interpretation of the results. The algorithm has the potential to be valuable for practitioners managing option positions on a single underlying asset.

Figure 4 . 1 :

 41 Figure 4.1: Optimal values of δ1 3 (right) and δ2 3 (left) as a function of δ 1 2 , δ 2 2 .

Figure 4 . 2 :

 42 Figure 4.2: Optimal values of δ1 3 (right) and δ2 3 (left) as a function of δ 1 2 , δ 2 2. Proportional cost for trading S 1 is 1%, while the proportional cost for trading S 2 has been raised to 10%.

Table 1 :

 1 The hedging strategies δ and δ.

	i 1 ) as the hedging

Table 2 :

 2 , have come to similar Comparison between hedging strategies δ and δ.

	Strategy δ	Strategy	δ

Table 3 :

 3 Parameters of the Heston model and of the put option.

	Meaning	Value
	S 0 Initial spot value	100
	r Risk free i.r.	0, 0.1
	µ Drift	0, 0.1
	v 0 Vol at beginning	0.04
	α Mean reversion speed	1
	b Long run vol	0.04
	σ Vol of Vol	2
	ρ BW Correlation	-0.7
	T Maturity	1.0
	K Strike	100
	Symbol Meaning		Value
	M 0 Number of MC paths for computing quantiles	10 6
	M Number of MC paths for backward step		10 3 , 10 4 , 10 5
	M 1 Number of MC paths for out-of-sample step	10 6
	N Number of trading dates		4, 8, 16
	N T Number of time-steps for simulation		16
	N S Number of points to discretyze S		15
	N I Number of points to discretyze the integral of V	2
	N δ 1 Number of points to discretyze δ 1 n-1		3
	N δ 2 Number of points to discretyze δ 2 n-1		3

Table 4 :

 4 Backward Hedging parameters.

	Symbol Meaning	Value
	M 0 Number of MC paths for computing quantiles	10 6
	M Number of MC paths for backward step	10 3 , 10 4 , 10 5
	N Number of trading dates	4, 8, 16
	N T Time-steps for simulation	64
	N N Number of neurons per layer	17
	N L Number of hidden layers	2
	σ(x) Activation function	max (x, 0)

Table 5 :

 5 Deep Hedging parameters.

Table 6 :

 6 Prices for a European put option with r = µ = 0, while minimizing loss based on mean squared error.

			ε = 0			ε = 1%
	Benchmark		3.7165			-
	M N	4	8	16	4	8	16
	10 3	3.7360	3.7513	3.8066	5.2907	5.6603	6.9647
		(8)	(16)	(49)	(86)	(382)	(1618)
	BH 10 4	3.7022	3.7172	3.6958	5.0861	5.3101	5.5282
		(12)	(35)	(111)	(407)	(2377)	(10195)
	10 5	3.7081	3.7038	3.7057	5.1669	5.1404	5.3706
		(42)	(197)	(872)	(3336)	(20906)	(88598)

Table 7 :

 7 Prices for an American put option with r = µ = 0, while minimizing loss based on mean squared error.

				ε = 0			ε = 1%
	Benchmark		2.1334			-
		M N	4	8	16	4	8	16
		10 3	2.1537	2.1470	2.2142	2.9850	3.1667	3.9847
			(6)	(10)	(18)	(96)	(391)	(1499)
	BH 10 4	2.1392	2.1547	2.1412	3.0313	3.0702	3.0806
			(8)	(14)	(31)	(364)	(2072)	(8895)
		10 5	2.1462	2.1401	2.1410	2.9695	2.8845	2.9150
			(16)	(43)	(147)	(2674)	(17137)	(74737)
		10 3	1.3124	1.3044	1.1886	2.3709	3.9069	5.2301
	DH	10 4	(208) 1.9725	(396) 1.9670	(677) 1.7516	(440) 3.4638	(506) 4.7569	(800) 4.6196
			(944)	(1981)	(4399)	(900)	(2427)	(4678)
		10 5	2.1270	2.0983	2.1339	3.2011	3.9313	4.1068
			(11319)	(20367)	(41625)	(27554)	(64813)	(85591)

Table 8 :

 8 Prices for a European put option with r = µ = 0.1, while minimizing loss based on mean squared error.

			ε = 0			ε = 1%
	Benchmark		2.8355			-
	M N	4	8	16	4	8	16
	10 3	2.6580	2.7427	2.8339	3.5394	3.8889	4.7330
		(6)	(10)	(20)	(87)	(387)	(1423)
	BH 10 4	2.6238	2.7503	2.7745	3.4659	3.8864	3.7912
		(8)	(16)	(39)	(444)	(2945)	(10535)
	10 5	2.6264	2.7376	2.7730	3.5296	3.5533	3.6406
		(20)	(64)	(203)	(4140)	(27944)	(93940)

Table 9 :

 9 Prices for an American put option with r = µ = 0.1, while minimizing loss based on mean squared error.

Table 10 :

 10 Hedging cost for a European put option with r = µ = 0, while minimizing loss based on CVaR.

			ε = 0			ε = 1%
	M N	4	8	16	4	8	16
	10 3	10.4560	8.7006	7.6962	11.5435	10.3831	9.3273
		(10)	(14)	(27)	(119)	(499)	(2087)
	BH 10 4	10.1111	8.6462	7.4082	11.5382	10.3588	9.3308
		(35)	(73)	(148)	(476)	(2733)	(11618)
	10 5	10.0944	8.6724	7.4855	11.5117	10.3395	9.3630
		(355)	(759)	(1625)	(3721)	(21980)	(91381)

Table 11 :

 11 Hedging cost for an American put option with r = µ = 0, while minimizing loss based on CVaR.

			ε = 0			ε = 1%
	M N	4	8	16	4	8	16
	10 3	6.5798	5.3969	4.4896	7.4402	6.1939	5.3032
		(13)	(13)	(23)	(139)	(497)	(1957)
	BH 10 4	6.5619	5.2388	4.3412	7.3739	6.0801	5.3044
		(35)	(70)	(140)	(448)	(2443)	(10104)
	10 5	6.5139	5.2283	4.3347	7.3460	6.0491	5.3076
		(359)	(803)	(1545)	(3354)	(20161)	(79368)
	10 3	18.1152	23.8998	24.1038	24.5469	27.4145	29.8575
		(388)	(522)	(1530)	(424)	(802)	(1689)
	DH 10 4	17.5067	19.2940	18.3218	18.2993	18.8594	21.2792
		(1620)	(1818)	(3710)	(2031)	(1892)	(5720)
	10 5	10.2705	12.2695	8.1678	10.3475	8.1222	7.9165
		(26602)	(11261)	(67796)	(42343)	(53501)	(105889)
	10 6	6.6167	5.8382	5.3762	6.9460	5.9089	5.2575
		(22974)	(31135)	(42442)	(19964)	(62982)	(110851)

Table 12 :

 12 Hedging cost for a European put option with r = µ = 0.1, while minimizing loss based on CVaR.

			ε = 0			ε = 1%
	M N	4	8	16	4	8	16
	10 3	8.2635	7.0993	6.0820	8.7475	7.7933	7.3570
		(9)	(14)	(30)	(117)	(438)	(1729)
	BH 10 4	8.3280	6.9922	6.0175	8.6478	7.9860	7.4585
		(34)	(71)	(143)	(454)	(2332)	(9742)
	10 5	8.2982	6.9864	5.9918	8.8270	8.0841	7.5112
		(341)	(731)	(1530)	(3539)	(19797)	(77681)

Algorithm 1 Backward Hedging for European options

Set up volatility tree T by following [START_REF] Goudenège | Gaussian process regression for pricing variable annuities with stochastic volatility and interest rate[END_REF] if ε > 0 then run this algorithm with ε = 0 to obtain the parameters for the 5D grids Gn for n = 1, . . . , N -1 else simulate M0 random paths S 1 n,m , In,m for k = 1, . . . , N for n = N -1, . . . , 0 do compute quantiles from S1 n,m , Ĩn,m and set the 3D grid Gn from quantiles and T end for for n = 0, . . . , N , the payoff φ S 1 T m , and the hedging error

e -rT and, for n = 1, . . . , N -1, the quantiles of S 1 tn , It n , δ1 n and δ2 n based on the out of sample simulations end if

Algorithm 2 Backward Hedging for American options

Set up volatility tree T by following [START_REF] Goudenège | Gaussian process regression for pricing variable annuities with stochastic volatility and interest rate[END_REF] if ε > 0 then run this algorithm with ε = 0 to obtain the parameters for the 5D grids Gn for n = 1, . . . , N -1 else simulate M0 random paths S 1 n,m , In,m for k = 1, . . . , N for n = N -1, . . . , 0 do compute quantiles from S 1 n,m , In,m and set the 3D grid Gn from quantiles and T end for end if for n = N -1, . . . , 0 do for all points G of Gn (if