RECOVERING PROBABILITY FUNCTIONS WITH FOURIER SERIES

ABSTRACT

The COS method was introduced in Fang & Oosterlee (2008) and then was applied to pricing a variety of stock options for continuous random variables. This paper adapts the Fourier-cosine series (COS) method to recover discrete probability mass functions. We approximate mixture and compound probability distributions with cosine series. Enormous precision and computational speed are the qualities of the function estimates here obtained. We also develop the pricing framework to trade derivatives subject to discrete random variables. We apply the method to calculate, for the first time, the price of an interest rate derivative of recent vintage introduced in the Brazilian financial market. Parameter calibration confirms the quality of the model.

Keywords:
probability functions; Fourier series; COS method

1 INTRODUCTION

There are situations in which the probability density or mass function is not available in closed form. Classical distributions (e.g. normal and lognormal densities or Poisson distribution) which lead to closed form expressions for a variety of problems, sometimes are not capable to reproduce the behavior seen in real-problem data. Fang & Oosterlee (2008) introduced a method intended to recover the probability density functions in a quasi-analytical way. Given the characteristic function of a continuously distributed random variable, the authors showed that it is possible - using the Euler’s identity, to approximate pointwise the probability density function via the Fourier-cosine series.

In several numerical experiments, these authors showed that the convergence rate of the COS method is exponential and its computational complexity linear. Oosterlee & Grzelak (2019) is a good text in which the COS method is explained in detail. Many works have recently extended the original COS model and applied it to solve a variety of option pricing problems under classical models. Zhang & Oosterlee (2013) showed how to use the COS method to pricing arithmetic and geometric Asian options under exponential Lévy processes. Zhang & Feng (2018) found the price of American put options under the double Heston (1993) model using the COS method. Zhang et al. (2012) analyzed the efficiency properties of pricing commodity options with early exercise. Ballotta et al. (2018) employed the COS method in multivariate structural models to pricing credit default using Lévy processes. Tour et al. (2018) applied the Fourier cosine expansion method to calculate the price of several options under regime-switching models. Have & Oosterlee (2018) used the COS method for option valuation under the Stochastic Alpha Beta Rho (SABR) dynamics. The pricing of forward starting options with stochastic volatility and jumps was considered in Zhang & Geng (2017). Crisóstomo (2018) compared the CPU effort of the COS method against six other Fourier-based schemes and concluded that it is notably the fastest. Newer methods have appeared in the literature aiming the pricing of European options. Notable examples are Ortiz-Gracia & Oosterlee (2013), Ortiz-Gracia & Oosterlee (2016) and Kumar et al. (2019), the first two based on wavelet techniques.

In spite of the vast existing literature concerning the above subject matter, to the best of our knowledge, except for da Silva et al. (2019), there is no adaptation of the COS method to interest rate derivatives pricing problems. In fact, there is a lack in the literature of a general fast and accurate pricing method for interest rate options in distinct environments. We introduce the scenario where jumps only occur at deterministic prescribed times with a restrict number of jump sizes. This addresses a very practical issue concerning the Monetary Policy Committee meetings and the prescribed dates when the committees change their basic interest rates.

We extend the method to recover probability mass functions for discrete random variables and show that the Fourier-cosine series with a finite number of terms is exact for some cases. In Section 2.1 we present some examples for mixture and compound probability distributions to show how the method can be useful when the probability density or mass function is not easy to manipulate, its cumulative function has no analytical solution and/or the function itself does not exist in an explicit form. In Section 3 we applied our results in the context of the interest rate financial markets. We develop the framewok to use the COS method to calculate the price of the Monetary Policy Committee (COPOM) Option using, for the first time in the financial literature, the (discrete) Skellam distribution, more precisely, a modified version of the Skellam probability distribution. Finally, we show that the model matches the market data.

2 FOURIER COSINE-SERIES FOR DISCRETE RANDOM VARIABLES

This is an extension of Fang & Oosterlee (2008) for discrete random variables. Consider a function f taking values in ℝ and defined on a discrete partition of [0, π], namely, $\Omega =\{0,1\delta ,...,n\delta =\pi \}$, with n a positive integer. Then f can be expressed by the finite Fourier-cosine series

$f\left(k\right)=\frac{a_0}{2}+\sum _{j=1}^n{a}_j\cos \left(jk\right),k\in \Omega$ (1)

where:

$a_j=\frac{2}{\pi }\sum _{k\in \Omega }f\left(k\right)\cos \left(jk\right),j=0,1,2,...,n.$ (2)

For functions supported in a partition of any arbitrary interval [a, b], namely, $\hat{\Omega }=\{a,a+\hat{\delta },...,a+n\hat{\delta }=b\},\hat{\delta }=\delta \frac{b-a}{\pi }$, a change of variable $k=\pi \frac{x-a}{b-a}$ is considered. Then, the Fourier-cosine series expansion of f in $\hat{\Omega }$ is

$f\left(x\right)=\frac{a_0}{2}+\sum _{j=1}^n{a}_j\cos \left(j\pi \frac{x-a}{b-a}\right),x\in \hat{\Omega }$ (3)

where:

$a_j=\frac{2}{b-a}\sum _{x\in \hat{\Omega }}f\left(x\right)\cos \left(j\pi \frac{x-a}{b-a}\right),j=0,1,2,...,n.$ (4)

By the Euler’s identity, the coefficients of the Fourier-cosine expansion of f are

$$\begin{gathered} a_j=\frac{2}{b-a}\sum _{x\in \hat{\Omega }}f\left(x\right)\mathfrak{R}\left[e^{\left(ij\pi \frac{x-a}{b-a}\right)}\right] \\ =\frac{2}{b-a}\mathfrak{R}\left[e^{\left(-ij\pi \frac{a}{b-a}\right)}\sum _{x\in \hat{\Omega }}f\left(x\right)e^{\left(ij\pi \frac{x}{b-a}\right)}\right]. \end{gathered}$$

Let X be a discrete random variable taking values in the extension of $\hat{\Omega }$, i.e., ${\hat{\Omega }}_e=\{...,a-\hat{\delta },a,a+\hat{\delta },...,a+n\hat{\delta }=b,b+\hat{\delta },...\}$. If f : $f:{\hat{\Omega }}_e\to \mathrm{ℝ}$ is the probability mass function of X, then

$a_j\approx \frac{2}{b-a}\mathfrak{R}\left[e^{\left(-ij\pi \frac{a}{b-a}\right)}\hat{f}\left(\frac{j\pi }{b-a}\right)\right]\triangleq A_j,$ (5)

where $\hat{f}$ is the characteristic function of X, that is

$\hat{f}\left(u\right)=\sum _{x\in {\hat{\Omega }}_e}{e}^{\left(ixu\right)}f\left(x\right),$ (6)

which is approximated in a more restricted domain $\hat{\Omega }$ by

$\sum _{x\in \hat{\Omega }}{e}^{\left(ixu\right)}f\left(x\right).$ (7)

Therefore, the approximation of f is given by the following Fourier-cosine series

$f\left(x\right)\approx \frac{A_0}{2}+\sum _{j=1}^{\hat{n}}{A}_j\cos \left(j\pi \frac{x-a}{b-a}\right),x\in \hat{\Omega },$ (8)

for $\hat{n}$ equal to the greatest integer less than $f\left(y\right)=0,\forall y\in \left\{{\hat{\Omega }}_e\setminus \hat{\Omega }\right\}$, then a _j = A _j and (8) is exact.

The choices of the limits a and b for the approximation (8) be good were proposed in Fang & Oosterlee (2008) as follows:

$a=c_1-L\sqrt{c_2+\sqrt{c_4}},b=c_1+L\sqrt{c_2+\sqrt{c_4}}$ (9)

with L = 10. The coefficients c _k used in (9) are the k-th cumulant function of X and they are given by

$c_k=\frac{1}{i^k}\frac{d^k}{d{u}^k}h\left(u\right){|}_{u=0}.$ (10)

The cumulant generating function h(u) that appears in (10) is given by

$h\left(u\right)=\ln \mathbb{E}\left[e^{\left(iuX\right)}\right].$ (11)

Remind that the domain of f, typically, is not [a, b]. This interval is chosen according to (9) in order to capture as much probability as possible from f. In Junike & Pankrashkin (2022) the authors showed how to obtain the truncation interval using Markov’s inequality ensuring convergence of the COS method within a predefined error tolerance. This could be an alternative to define a and b. Notice however that, depending on the domain of the function, we can reach the exact result in some few steps of the convergence series, preserving the approach according to Fang & Oosterlee (2008) as we do in (9).

The discrete case introduced above was not addressed in Fang & Oosterlee (2008). The (huge) paper Abate & Whitt (1992) gives a glimpse in this subject matter. However, a discussion of discrete cumulative distribution function appears in Ruijter et al. (2015). The authors use the continuous version of the COS method to calculate the cumulative probability, where the Gibbs phenomenon surges. They develop spectral filter to ensure convergence. In our discrete version of the COS method, we have

$\mathrm{\mathbb{P}}(K\le k)=\sum _{t=k_{min}}^k\mathrm{\mathbb{P}}(K=t).$ (12)

Example 1. (Binomial distribution) The probability of getting k successes in n independent Bernoulli trials is given by the binomial distribution with parameters (n, p), namely

$\mathrm{\mathbb{P}}(X=k)=\left(\genfrac{}{}{0ex}{}{n}{k}\right)p^k(1-p{)}^{n-k},$ (13)

where

$\left(\genfrac{}{}{0ex}{}{n}{k}\right)=\frac{n!}{k!(n-k)!}.$

Its characteristic function is given by

$\hat{f}(u;n,p)={\left(p{\mathrm{e}}^{\mathrm{i}u}+(1-p)\right)}^n.$ (14)

In order to calculate binomial probabilities with the COS method, i.e. use (8) as an approximation for (13), we need to choose the domain [a, b]. According to (11), we have the cumulant generating function of the binomial distribution given by

$h\left(u\right)=n\ln \left(p{\mathrm{e}}^{\mathrm{i}u}+(1-p)\right),$ (15)

and the cumulants

$c_1=np,c_2=\frac{np(1-p)}{{\left(p+(1-p)\right)}^2},c_4=\frac{np(1-p)\left(p^2-4p(1-p)+(1-p{)}^2\right)}{{\left(p+(1-p)\right)}^4}.$

InFigure 1we show the approximation of the binomial probability function with parameters p = 0.05 and n = 80 via the Fourier Series.

Figure 1
Binomial probability mass function and cumulative distribution.

2.1 Fourier-cosine series for mixture and compound probability distributions

The approximation given by (8) is particularly useful when the probability mass function is not easy to manipulate, its cumulative function has no analytical solution and/or the function itself does not exist in an explicit form. In such cases, we may resort to the Fourier cosine series given that we know the corresponding characteristic function.

In order to illustratethis, we calculate the mass functions of two probability distributions with non existing analytical formulas, but with existing analytical characteristic functions, namely, the Neyman type A distribution and the Skellam distribution, both represented by infinite sums.

A factorial function, which is difficult to obtain for large numbers, is also present as will be seen below. Moreover, for the latter, a gamma function must be solved for each term.

Example 2.(Neyman type A distribution) Consider a Poisson distribution for which the Poisson parameter Θ = ϕθ, where ϕ is a constant and θ has a Poisson distribution with parameter λ. The outcome distribution is known as Neyman type A probability distribution. This probability distribution found applications in biological systems and in the natural disaster modeling ((Martin & Katti, 1962) and (Ozel & Turkan, 2022)). Detailed discussion about mixture distributions can be found inJohnson et al. (2006).

The characteristic function of the Neyman type A probability distribution is given by

$\hat{f}(u;\lambda ,\varphi )=e^{\left[\lambda \left(e^{\varphi \left(e^{iu}-1\right)}-1\right)\right]}.$ (16)

According to (10), we have

$h\left(u\right)=\left[\lambda \left(e^{\varphi \left(e^{iu}-1\right)}-1\right)\right],$ (17)

and

$c_1=\lambda \varphi c_2=\lambda \varphi (1+\varphi )c_4=\lambda \varphi (1+7\varphi +6{\varphi }^2+{\varphi }^3).$

InFigure 2we show the approximation of the Neyman type A probability function via the Fourier Series. We highlight the multimodality of the distribution. The local peaks in 0, 25, 50, 76, 103, 129, 153 and 173 for the parameters λ = 25 and ϕ = 25 was discussed inShenton & Bowman (1967).

Figure 2
Neyman type A probability mass function.

An advantage of the solution via Fourier series is that the distribution is well approximated with a finite number of cosine terms in contrast to the Neyman’s type A analytical probability mass function given by

$\mathrm{\mathbb{P}}\{X=x\}=\frac{e^{-\lambda }{\varphi }^x}{x!}\sum _{j=0}^{\infty }\frac{(\lambda e^{-\varphi }{)}^j{j}^x}{j!},x=0,1,...,$ (18)

which has an infinite sum for each x and a factorial function. The latter is a computational challenge for large values of x.

Example 3.(Skellam distribution) The Skellam distribution is the discrete probability distribution of the difference N₁ − N ₂ of two statistically independent random variables N ₁ and N ₂ , each Poisson-distributed with respective expected values µ ₁ and µ ₂ and support on {..., −2, −1, 0, 1, 2, ...}.

In the study of cameras, the Skellam distribution is used to measure the intensity difference of pixels in the spatial and temporal domain as the photons are Poisson distributed ( Hwang et al., 2007 ). Application of the Skellam distribution is found in the investigation of football outcomes with modelling the number of goals scored by each team ( Karlis & Ntzoufras, 2009 ).

The characteristic function of the Skellam distribution is

$\hat{f}(u;{\mu }_1,{\mu }_2)=e^{-({\mu }_1+{\mu }_2)+{\mu }_1{e}^{iu}+{\mu }_2{e}^{-iu}}$ (19)

According to (10), we have

$h\left(u\right)=-({\mu }_1+{\mu }_2)+{\mu }_1{e}^{iu}+{\mu }_2{e}^{-iu},$ (20)

with

$c_1={\mu }_1-{\mu }_2,c_2={\mu }_1+{\mu }_2+({\mu }_1-{\mu }_2{)}^2$

and

$$\begin{gathered} c_4={\mu }_1^4+6{\mu }_1^3-{\mu }_1^2\left(4{\mu }_1{\mu }_2-7\right)-{\mu }_1\left(6{\mu }_1{\mu }_2-1\right)+2{\mu }_1{\mu }_2\left(3{\mu }_1{\mu }_2-1\right) \\ -{\mu }_2\left(6{\mu }_1{\mu }_2-1\right)-{\mu }_2^2\left(4{\mu }_1{\mu }_2-7\right)+6{\mu }_2^3+{\mu }_2^4. \end{gathered}$$

The probability mass function for the Skellam distribution for a difference $K=N_1-N_2$ is given by

$p(k;{\mu }_1,{\mu }_2)=\mathrm{\mathbb{P}}\{K=k\}=e^{-({\mu }_1+{\mu }_2)}{\left(\frac{{\mu }_1}{{\mu }_2}\right)}^{k/2}{I}_k\left(2\sqrt{{\mu }_1{\mu }_2}\right).$ (21)

where I_k (·) is the modified Bessel function of first kind, which is an ordinary differential equation with the following solution:

$I_k\left(z\right)={\left(\frac{z}{2}\right)}^k\sum _{i=0}^{\infty }\frac{{\left(\frac{z^2}{4}\right)}^i}{i!\Gamma (k+i+1)}.$ (22)

where

$\Gamma \left(x\right)=\underset{0}{\overset{\infty }{\int }}{s}^{x-1}{e}^{-s}ds.$ (23)

An advantage of the solution via Fourier series is that the distribution is well approximated with a small finite number of cosine terms in contrast to the analytical probability mass function given by (21), in which case a computational problem arises for large values of k in the modified Bessel function (23). Using MATLAB R2022a, the COS method is 4 times faster than the implementation of (21) for calculating the Skellam probability with the same number of terms. InFigure 3we show the approximation of the Skellam probability function via the Fourier Series with parameters µ ₁ = 25 and µ ₂ = 5. InFigure 4we show the convergence of the COS method when the number of summation terms increases and inFigure 5it is exhibit the error of the approximations. It is noteworthy that the COS method converges earlier than the analytical solution.

Figure 3
Skellam probability mass function.

Figure 4
Skellam distribution probability varying the number of summation terms.

Figure 5
Absolute error varying the number of summation terms.

Note that,$\mathbb{E}\left[K\right]={\mu }_1-{\mu }_2$. Then, it is computationally complex to calculate ℙ{K = k} with (21) for large values of 𝔼[K]. On the other hand, the level of parameters µ ₁ and µ ₂ has no impact in solving the probability via Fourier series.

3 PRICING OPTIONS UNDERLINED BY DISCRETE RANDOM VARIABLES

Financial instruments are usually modeled through continuous random variables. Black & Scholes (1973) won a nobel prize for showing that the price of a contingent claim is the expected value of the discounted payoff given the present information and that it results in a riskless portfolio. They used the geometric Brownian motion as the underlying asset model. Vasicek (1977) used a Gaussian structure to calculate bond prices. Duffie & Singleton (2003) introduced the con cept of affine-jump diffusion models to the context of credit risk, which was extended to other areas of finance.

Continuous stochastic processes dominate all areas of finance and discrete methods to solving time continuous problems arise when no analytical solution exists to such continuous problems. The research thrust concerning a discrete approach to time continuous problems can be perceived, for instance, via Duffy (2006) and Wilmott (2006) which present various finite difference methods used to find the price of financial derivatives. A particular finite difference method developed to calculate the price of interest rate derivatives can be found in da Silva et al. (2016). Monte Carlo simulations are shown in Glasserman (2004) and binomial trees are presented in Hull (2009).

Numerical methods such as finite differences, binomial trees and other approximations of continuous variables create discrete approximations that, in the limit of small step-sizes, converge to the continuous solutions of these variables. On the other side, inherently discrete variables, like jump processes and discrete distributions will have a discrete reference distribution.

Now, scant attention is given for problems that are inherently of discrete time nature. With this in mind, our aim here is threefold:

• we show how to calculate the price of a contingent claim where the underlying asset is inherently a discrete random variable;

• we introduce a new probability distribution to deal with the possible changes in the basic interest rate due to the Monetary Policy Committee meetings;

• we calibrate the interest rate discrete jump model to market probabilities.

3.1 Pricing discrete options

Let x(t) be the observed discrete random variable at time t, with values in ${\hat{\Omega }}_e$, representing the underlying source of risk of an European call option maturing at T. In the financial jargon, European options can only be exercised at maturity (Hull, 2009). Denote by f the probability mass function of x(T ), defined in ${\hat{\Omega }}_e$. Thus, the price of this option at time t is

$C(t,T)=E\left[g\right(x\left(T\right)\left)x\right(t\left)\right]=\sum _{x\in {\hat{\Omega }}_e}g\left(x\right)f\left(x\right|x\left(t\right))dx,$ (24)

where g(x(T)) is the discounted payoff function of the option. Truncating f in the interval [a, b] we have

$C(t,T)\approx \sum _{x\in \hat{\Omega }}g\left(x\right)f\left(x\right|x\left(t\right)).$ (25)

Therefore, using f (x) as in (8) we have

$C(t,T)\approx \frac{A_0}{2}\sum _{x\in \hat{\Omega }}g\left(x\right)+\sum _{j=1}^{\hat{n}}{A}_j\sum _{x\in \hat{\Omega }}g\left(x\right)\cos \left(j\pi \frac{x-a}{b-a}\right).$ (26)

Hence, the series approximation of the option price is given by

$C(t,T)\approx \frac{A_0{B}_0}{2}+\sum _{j=1}^{\overline{n}}{A}_j{B}_j$ (27)

where the A _j coefficients are given by (5) and

$B_j=\sum _{x\in \hat{\Omega }}g\left(x\right)\cos \left(j\pi \frac{x-a}{b-a}\right).$ (28)

3.2 Pricing COPOM options

The interest rate derivatives market differs across countries. The existing singularities in the Brazilian market stirred the development of particular derivatives contracts. The most important is the IDI option discussed in da Silva et al. (2016) and in Carreira & Brostowicz (2016). Among them, we have a singular kind of interest rate option, namely the Monetary Policy Committee (COPOM) Option, a derivative of recent vintage introduced in the market. COPOM option is a cash-or-nothing product which pays a fixed amount if the discrete movement of the Selic Target Rate X defined in a given Copom meeting at T is equal to the traded change K and pays nothing otherwise. To the best of our knowledge, there is no stochastic formula to calculate the price of the COPOM option in the existing literature. The Copom option price is then given by

$C(t,T)={\mathrm{e}}^{-r\left(t\right)(T-t)}\mathbb{E}\left[g\left(X\right(T\left)\right)\right],$ (29)

where r(t) is the actual, whence observed, interest rate and

$g\left(X\right(T\left)\right)=1_{\left\{X\right(T)=K\}}.$

Then

$C(t,T)={\mathrm{e}}^{-r\left(t\right)(T-t)}\mathrm{\mathbb{P}}\left[X\left(T\right)=K\right],$ (30)

where $\mathrm{\mathbb{P}}\left[X\left(T\right)=K\right]$ is the probability that the movement of the Selic Target Rate equals the strike in the meeting occurring at T.

The monetary authority holds eight regularly scheduled meetings during the year. The scheduled central bank announcements move the benchmark rate discretely. Figure 6 shows the path of the basic target interest rate of the last two decades.

Figure 6
Brazilian SELIC rate.

In order to calculate the price of the Copom option, we need to calculate the probability of a discrete random variable which has support on the set of multiples of 0.25% and discount it to the present value.

In the theorem below, we introduce a modified version of the Skellam probability distribution shown in Section 2.1 to characterize the jump sizes of the SELIC target interest rate. The domain of the original Skellam distribution is the integer numbers {..., −2, −1, 0, 1, 2, ...}, as it is constructed as a difference of two independent Poisson processes.

Theorem 1. Let X be a discrete random variable governing by the Skellam distribution characterized by the following moment generating function:

$M_X\left(t\right)=e^{-({\mu }_1+{\mu }_2)+{\mu }_1{e}^t+{\mu }_2{e}^{-t}}$ (31)

Then the moment generating function of the modified Skellam distribution of $Y=aX+b$ is given by

$M_Y\left(t\right)=e^{-({\mu }_1+{\mu }_2)+bt+{\mu }_1{e}^{at}+{\mu }_2{e}^{-at}}.$ (32)

Proof. The moment generating function of Y is

$M_Y\left(t\right)=\mathbb{E}\left[e^{Yt}\right]=\mathbb{E}\left[e^{aXt+bt}\right]=e^{bt}\mathbb{E}\left[e^{aXt}\right]=e^{bt}{M}_X\left(at\right).$

So, it follows that

$M_Y\left(t\right)=e^{-({\mu }_1+{\mu }_2)+bt+{\mu }_1{e}^{at}+{\mu }_2{e}^{-at}}$ (33)

is the moment generating function of the modified Skellam probability distribution.

□

The modified Skellam distribution refers to the distribution of $Y=aX+b$, where X has a Skellam distribution. Clearly the domain of Y now is

$\left\{...,-2a+b,-a+b,b,a+b,2a+b,...\right\}.$

So, Y is not of Skellam type any more, except for (a, b) = (1, 0).

Corollary 1.Let a$a=\frac{1}{400}$and b = 0. Then, the modified Skellam distribution has the domain equal to the jumps of the SELIC target rates. The moment generating function is given by

$M_Y\left(t\right)=e^{-({\mu }_1+{\mu }_2)+{\mu }_1{e}^{\frac{t}{400}}+{\mu }_2{e}^{-\frac{t}{400}}}$ (34)

Proof. It is easy to check that $a=\frac{1}{400}$ and b = 0 in (33) implies (34). The domain is given by

$\{...,-0.50\%,-0.25\%,0,0.25\%,0.50\%,...\}.$

□

To the best of our knowledge, there is no other discrete random variable with the above properties. The example below shows that the modified Skellam distribution is properly applied in the pricing procedure.

Example 4.Suppose that the moment generating function parameters of the modified Skellam probability distribution (34) is µ ₁ = 4 and µ ₂ = 0.5. Let T − t = 0.1 be the maturity of the COPOM option, r(t) = 0.02 be the actual interest rate and K = 0.75% be the strike. So, the price of the COPOM option is given by (30)

$C(t,T)=f(0.0075)={\mathrm{e}}^{-(0.02)0.1}\left[\frac{A_0}{2}+\sum _{j=1}^n{A}_j\cos \left(j\pi \frac{0.0075-a}{b-a}\right)\right]=0.1906,$ (35)

i.e. the discounted probability of a change of +0.75% in the target interest rate given by the modified Skellam distribution at the next Copom meeting is 19.06%. Note that$f\left(y\right)=0,\forall y\in \left\{{\hat{\Omega }}_e\setminus \hat{\Omega }\right\}$, then a _j = A _j and the price is exact. Then, the price of the Copom option is given by the discounted probability given by (35) times the size of the contract.

The coefficients A _j are given by (5), where the characteristic function is M(t) in (34), namely

$\hat{f}(t;{\mu }_1,{\mu }_2)={\mathrm{e}}^{4{\mathrm{e}}^{\frac{\mathrm{i}t}{400}}+0,5{\mathrm{e}}^{-\frac{\mathrm{i}t}{400}}-4,5}.$ (36)

According to (11), the cumulant generating function is

$h\left(t\right)=4{\mathrm{e}}^{\frac{\mathrm{i}t}{400}}+0,5{\mathrm{e}}^{-\frac{\mathrm{i}t}{400}}-4,5.$ (37)

So, the cumulants are

$c_1=\frac{3,5}{400},c_2=\frac{4,5}{{400}^2},c_4=\frac{4,5}{{400}^4}.$ (38)

Actually, the Selic target interest rate X(T) derives from the short rate r(t) which, in turn, is governed by a stochastic differential equation (SDE). Now, in Genaro & Avellaneda (2018) the authors advocate that, in an inflation economy, interest rate changes are decomposed into fluctuations of the short-term interest rate between Central Bank meetings and deterministic timed jumps following the meetings. Hence, to dully model the short interest rate r(t), a time deterministic jump term should enter the SDE. Namely, we shaped this term as follows.

$\sum _{j=1}^p{1}_{\{t={\tau }_j\}}{J}_{{\tau }_j}dN\left(t\right),$ (39)

where N(t) is a counting process, the summation carries out the activation of the jumps at prescribed deterministic times τ _j , j = 1, 2, ..., and $J_{{\tau }_j}$ are random variables that represent the jump amplitudes of the short term interest rate target r(t).

In order to accomplish this realistically, $J_{{\tau }_j}$ should take values in a prescribed discrete state space and, ultimately, be flexible enough in order to track the movement of the Selic target Rate X(T). The modified Skellam mass probability function is indeed a rare finding aimed to model the SDE jump term in the interest rate scenario, as it was confirmed in da Silva (2021). It is worth noticing that in the works of Genaro & Avellaneda (2018) jump intensities occur in the same prescribed state space as ours, however, our approach, via the Skellam distribution, encompasses any affine model, while theirs are restricted to the (affine) Black model (Black (1976)).

3.3 Calibration

Calibration of interest rate models can be seen as an optimization problem, where the parameters of the pricing function are determined to fit the observable prices.

The above modified Skellam probability moment generating function (34) have two parameters µ ₁ and µ ₂. In order to preclude arbitrage opportunities, we have to choose the model parameters so that the theoretical prices match the prices of traded contracts. Kienitz & Wetterau (2012) provide a deep review of optimization techniques for parameter calibration.

We access the B3 market data to get the target interest rate jump probabilities for the subsequent COPOM meeting. The results are the following:

Table 1
COPOM option market data (08/30/2022).

Movement	Probability
0.00%	0.729
+0.25%	0.225
+0.50%	0.040
+0.75%	0.006

• Source: https://www.b3.com.br/pt br/market-data-e-indices/servicos-de-dados/datawise/dashboard-publico/opcao-de-copom.htm

We use the function fgoalattain of the MATLAB R2022a to search the optimal parameters of the modified Skellam probability distribution. The objective is to find the minimum distance between probabilities given by the modified Skellam probability distribution and the market probabilities shown in the table 1 subject to µ ₁ > 0 and µ ₂ > 0. The resulting probabilities and the corresponding discrepancies are shown in figures 7a and 7b, respectively.

Figure 7
Calibration results.

Matlab takes 0.085711 seconds to find µ ₁ = 0.3109 and µ ₂ = 0.0031. As can be seen, the error is of the order of 10⁻³. The computer used for all experiments has an Intel Core i5 CPU, 2.53GHz.

4 CONCLUSIONS AND WORKS TO COME

We extended the COS method to recover probability distributions of discrete random variables. Or else, assuming that the analytical solution of the random variable’s characteristic function exists, the method allows us to recover the corresponding probability mass function through a finite series. We show how precisely the cosine series, with few terms even, approximates the discrete, compounded and mixture probability distributions.

Concerning the financial niche, we developed the framework in order to calculate - via a finite sum - the price of a novel derivative product traded at the Brazilian market, the COPOM option. Contributing to this, we introduced for the first time in the literature the Skellam distribution and used a dully modified version of it to enter the stochastic differential equation. Such version ultimately lead us to a good representation of the target interest rate jump probabilities for the subsequent COPOM meeting.

We also show that the parameter calibration is effective and fast and that the theoretical probability distribution matches the one observed in the market.

Relying on Subsection 3.2, we have that the stochastic differential equations for the short term interest rate can be enhanced with jumps shaped with the modified Skellam distribution enabling more reliable prices.

Other interest rate derivatives can be investigated using our framework since prices of IDI options and Future DI options are also subjected to the movements of the target rate.

Acknowledgment

This work was supported by CNPq under the grant 140840/2015-0.

References

• ABATE J & WHITT W. 1992. The Fourier-series method for inverting transforms of probability distributions. Queueing Systems, pp. 5-88.
• BALLOTTA L, FUSAI G & MARAZZINA D. 2018. Integrated structural approach to credit value adjustment. European Journal of Operational Research, 272(3): 1143-1157.
• BLACK F. 1976. The pricing of commodity contracts. Journal of Financial Economics, 3(1-2): 167-179.
• BLACK F & SCHOLES M. 1973. The pricing of options and corporate liabilities. Journal of Political Economy, 81: 637-54.
• CARREIRA M & BROSTOWICZ R. 2016. Brazilian Derivatives and Securities: Pricing and Risk Management of FX and Interest-Rate Portfolios for Local and Global Markets. Palgrave Macmillan UK.
• CRISÓSTOMO R. 2018. Speed and biases of Fourier-based pricing choices: a numerical analysis. International Journal of Computer Mathematics, 95(8): 1565-1582.
• DA SILVA AJ. 2021. Fast Pricing of Path-Dependent Interest Rate Options with Jumps in Continuous and Discrete Time and Stochastic Volatility. Ph.D. thesis. National Laboratory for Scientific Computing - LNCC.
• DA SILVA AJ, BACZYNSKI J & BRAGANÇA JFS. 2019. Path-Dependent Interest Rate Option Pricing with Jumps and Stochastic Intensities. Lecture Notes in Computer Science, 11540: 710-716.
• DA SILVA AJ, BACZYNSKI J & VICENTE JVM. 2016. A new finite method for pricing and hedging fixed income derivatives: Comparative analysis and the case of an Asian option. Journal of Computational and Applied Mathematics, 297: 98-116.
• DUFFIE D & SINGLETON KJ. 2003. Credit risk: pricing, measurement, and management. Princeton University Press.
• DUFFY DJ. 2006. Finite Difference Methods in Financial Engineering: A Partial Differential Equation Approach. John Wiley & Sons.
• FANG F & OOSTERLEE CW. 2008. A Novel Pricing Method for European Options Based on Fourier-Cosine Series Expansions. SIAM Journal on Scientific Computing, 31(2): 826-848.
• GENARO AD & AVELLANEDA M. 2018. Pricing Interest Rate Derivatives Under Monetary Policy Changes. International Journal of Theoretical and Applied Finance, 21(6): 1850037.
• GLASSERMAN P. 2004. Monte Carlo Methods in Financial Engineering. Applications of mathematics : stochastic modelling and applied probability. Springer.
• HAVE Z & OOSTERLEE CW. 2018. The COS method for option valuation under the SABR dynamics. International Journal of Computer Mathematics , 95(2): 444-464.
• HESTON SL. 1993. A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options. The Review of Financial Studies, 6(2): 327-343.
• HULL JC. 2009. Options, Futures and Others Derivatives. 7th ed.. Pearson Prentice Hall.
• HWANG Y, KIM JS & KWEON IS. 2007. Sensor noise modeling using the Skellam distribution: Application to the color edge detection. In: 2007 IEEE Conference on Computer Vision and Pattern Recognition. pp. 1-8.
• JOHNSON NL, KEMP AW & KOTZ S. 2006. Univariate Discrete Distributions. 3rd ed.. Wiley Series in Probability and Statistics. Wiley.
• JUNIKE G & PANKRASHKIN K. 2022. Precise option pricing by the COS method-How to choose the truncation range. Applied Mathematics and Computation, 421: 126935.
• KARLIS D & NTZOUFRAS I. 2009. Bayesian modelling of football outcomes: using the Skellam’s distribution for the goal difference. IMA Journal of Management Mathematics, 20(2): 133-145.
• KIENITZ J & WETTERAU D. 2012. Financial Modelling: Theory, Implementation and Practice with MATLAB Source. The Wiley Finance Series. Wiley.
• KUMAR P, RODRIGO C, GASPAR FJ & OOSTERLEE CW. 2019. On Local Fourier Analysis of Multigrid Methods for PDEs with Jumping and Random Coefficients. SIAM Journal on Scientific Computing , 41(3): A1385-A1413.
• MARTIN DC & KATTI SK. 1962. Approximations to the Neyman Type A Distribution for Practical Problems. Biometrics, 18(3): 354-364. Available at: http://www.jstor.org/stable/2527477
  » http://www.jstor.org/stable/2527477
• OOSTERLEE CW & GRZELAK LA. 2019. Mathematical Modeling and Computation in Finance. World Scientific.
• ORTIZ-GRACIA L & OOSTERLEE CW. 2013. Robust Pricing of European Options with Wavelets and the Characteristic Function. SIAM Journal on Scientific Computing , 35(5): B1055-B1084.
• ORTIZ-GRACIA L & OOSTERLEE CW. 2016. A Highly Efficient Shannon Wavelet Inverse Fourier Technique for Pricing European Options. SIAM Journal on Scientific Computing , 38(1): B118-B143.
• OZEL G & TURKAN S. 2022. Neyman type A Distribution for the Natural Disasters and Related Casualties in Turkey. Journal of Data Science, 13(3): 533-550.
• RUIJTER MJ, VERSTEEGH M & OOSTERLEE C. 2015. On the application of spectral filters in a Fourier option pricing technique. Journal of Computational Finance, 19(1): 75-106.
• SHENTON LR & BOWMAN KO. 1967. Remarks on Large Sample Estimators for Some Discrete Distributions, Technometrics. Technometrics, 9(4): 587-598.
• TOUR G, THAKOOR N, KHALIQ AQM & TANGMAN DY. 2018. COS method for option pricing under a regime-switching model with time-changed Lévy processes. Quantitative Finance, 18(4): 673-692.
• VASICEK O. 1977. An Equilibrium Characterization of the Term Structure. Journal of Financial Economics , 5: 177-188.
• WILMOTT P. 2006. Paul Wilmott on Quantitative Finance . 2th ed.. Chichester: John Wiley & Sons.
• ZHANG B, GRZELAK LA & OOSTERLEE CW. 2012. Efficient pricing of commodity options with early-exercise under the Ornstein-Uhlenbeck process. Applied Numerical Mathematics, 62(2): 91-111.
• ZHANG B & OOSTERLEE CW. 2013. Efficient Pricing of European-Style Asian Options under Exponential Lévy Processes Based on Fourier Cosine Expansions. SIAM Journal on Financial Mathematics, 4(1): 399-426.
• ZHANG S & GENG J. 2017. Fourier-cosine method for pricing forward starting options with stochastic volatility and jumps. Communications in Statistics - Theory and Methods, 46(20): 9995-10004.
• ZHANG SM & FENG Y. 2018. American option pricing under the double Heston model based on asymptotic expansion. Quantitative Finance , pp. 1-16.

Publication Dates

• Publication in this collection
  13 Mar 2023
• Date of issue
  2023

History

• Received
  15 Sept 2022
• Accepted
  26 Dec 2022