The Java applet for pricing Asian options under Heston’s model using the new Ninomiya weak approximation scheme and quasi-Monte Carlo

The Java applet for pricing Asian options under Heston’s model using the new

Ninomiya weak approximation scheme and quasi-Monte Carlo

by

Boyko Vasilev

Kandidatarbete i matematik / tillämpad matematik

MÄLARDALEN UNIVERSITY UKK / TM

Box 883

Date: 2008-05-14 Project name:

The Java applet for pricing Asian options under Heston’s model using the new Ninomiya weak approximation scheme and quasi-Monte Carlo

Author(s):

Boyko Vasilev (boyko.vasilev@gmail.com) Supervisor(s): Anatoliy Malyarenko Examiner: Dimitri Silvestrov Comprising: 15 ECTS credits (hp)

This study is based on a new weak-approximation scheme for stochastic differential equa-tions applied to the Heston stochastic volatility model. The scheme was published by Ni-nomiya and NiNi-nomiya (2008) and is an extension of Kusuoka’s approximation scheme.

Ninomiya’s algorithm decomposes Kusuoka’s stochastic model into a set of ordinary dif-ferential equations with random coefficients and suggests several numerical optimisations for faster calculation.

The subject of this paper is a Java applet which calculates the price of an Asian option under the Heston model.

I would like to thank my supervisor, Anatoliy Malyarenko, for his guidance and help during my work on this thesis. I would also like to thank Michael Meyer from the Martingale project for his implementation of quasi-random generator.

The subject of this thesis work is building a Java Applet for calculating the price of Asian average strike call options. The price is calculated based on the Heston stochastic volatility model using a new weak approximation scheme by Ninomiya and Ninomiya (2008). The applet implements a quasi-Monte Carlo simulation method for faster convergence of the approximation algorithm. The user can specify the number of simulated paths in order to achieve the desired accuracy. The applet outputs the price of the option and a graphical representation of the change in price depending on the number of simulations. The applet was built with the Eclipse development environment and the Java 6 SDK. The free java library JFreeChart was used for plotting graphical output.

1 The Problem 6

1.1 The Black and Scholes Model . . . 6

1.2 Asian options . . . 6

2 Theoretical framework 8 2.1 Notation . . . 8

2.2 The Heston model . . . 8

2.3 The new Ninomiya weak approximation scheme . . . 9

3 Implementation of the algorithm 13 3.1 Monte Carlo integration . . . 13

3.2 quasi-Monte Carlo integration . . . 14

3.2.1 Low-discrepancy sequences . . . 15

3.2.2 Normally distributed quasi-random sequences . . . 16

3.3 The Java applet . . . 18

3.3.1 Input . . . 18

3.3.2 Output . . . 19

4 Analysis 21 4.1 Experiments . . . 21

4.2 Convergence and integration error . . . 25

4.3 Accuracy of the algorithm . . . 25

5 Conclusion 26

A Source-code listing of Hestonwa.java 28

B Source-code listing for Derivs.java 46

3.1 Integration of f(x)over[a, b] . . . 13

3.2 Scatter-plot of pseudo- and quasi-random sequences, N=100 . . . 15

3.3 Inverted uniform and native standard normal sequences, N =100 . . . 17

3.4 The applet interface . . . 18

3.5 Check for correctness of the correlation coefficient . . . 19

3.6 Check for correct inputs . . . 19

3.7 Check for the inequality 2αθ−β2 >0 . . . 19

4.1 Asian call option value depending on stock price, given K = 1.05, r = 0.05, σ =0.09, T =1, α=2, β =0.1, θ =0.09, ρ =0, time-steps = 48 and N =3000 simulations. . . 21

4.2 Asian call option value depending on strike price, given S = 1, r = 0.05, σ =0.09, T =1, α=2, β =0.1, θ =0.09, ρ =0, time-steps = 48 and N =3000 simulations. . . 22

4.3 Asian call option value depending on risk-free interest, given S =1, K =1.05, σ =0.09, T =1, α=2, β =0.1, θ =0.09, ρ =0, time-steps = 48 and N =3000 simulations. . . 22

4.4 Asian call option value depending on annualized volatility, given S =1, K = 1.05, r = 0.05, T = 1, α = 2, β = 0.1, θ = 0.09, ρ = 0, time-steps = 48 and N =3000 simulations. . . 22

4.5 Asian call option value depending on time to maturity, given S=1, K =1.05, r = 0.05, σ = 0.09, α = 2, β = 0.1, θ = 0.09, ρ = 0, time-steps = 48 and N =3000 simulations. . . 23

4.6 Asian call option value depending on speed of mean-reversion of the volatil-ity, given S = 1, K =1.05, r =0.05, σ =0.09, T =1, β =0.1, θ =0.09, ρ =0, time-steps = 48 and N =3000 simulations. . . 23

4.7 Asian call option value depending on the long-run mean value for volatility, given S = 1, K = 1.05, r = 0.05, σ = 0.09, T = 1, α = 2, β = 0.1, ρ = 0, time-steps = 48 and N =3000 simulations. . . 23

4.8 Asian call option value depending on the volatility of volatility, given S = 1, K = 1.05, r = 0.05, σ = 0.09, T = 1, α = 2, θ = 0.09, ρ = 0, time-steps = 48 and N =3000 simulations. . . 24

4.9 Asian call option value depending on correlation between stock price and volatility, given S = 1, K = 1.05, r = 0.05, σ = 0.09, T = 1, α = 2, β = 0.1, θ =0.09, time-steps = 48 and N =3000 simulations. . . 24

4.10 Asian call option value depending on number of time-steps, given S = 1, K = 1.05, r = 0.05, σ = 0.09, T = 1, α = 2, β = 0.1, θ = 0.09, ρ = 0 and N =3000 simulations. . . 24

Chapter 1

The Problem

1.1

The Black and Scholes Model

The highly-acknowledged Black and Scholes model was created in 1973 by Fischer Black and Myron Scholes. The model can price non-dividend paying options of European type and has been used widely for pricing securities since the 1970s. Despite its robustness and simple calculation, the model has proved empirically wrong. It does not fit the real market prices perfectly and gives only a rough approximation.

The reason for that is the way the model treats the volatility of the underlying asset. Even though the model assumes that the underlying’s price follows a geometrical Brownian motion (GBM), it uses implied volatility throughout the life of the security. This gives a good approximation for the value of at-the-money options, but is far from fair for in- and out-of-the money options. Careful examination of market data shows that Black and Scholes’ implied volatility for different strikes and maturities are not constant. Plotting the volatility surface from real market data shows that it forms a so called volatility smile, which explains why using one single implied volatility is not reliable. Karoui et al. (1998) shows that the original Black and Scholes formula is very inaccurate for pricing path-dependent stochastic volatility options.

1.2

Asian options

Asian options(also called average options) are becoming increasingly popular in times when instability reigns both the financial and commodity markets. Asian options are exotic op-tions and their underlying variable is not the asset price at maturity, but the average price over a period of time. This type of option has lower volatility and is therefore cheaper than an European option. Average options are attractive to buyers of commodities and energy goods because these good are usually traded at average prices over a period of time. Aver-age options also offer some protection against exchange-rate risks for exporters and multi-national corporations. Average options are less vulnerable to price manipulations, because it is more difficult to manipulate the average price of the underlying, than to affect it only at the point of maturity. Average options were traded for the first time in 1987 in Tokyo, where Bankers Trust developed a pricing formula for options based on the average price of crude oil. Because it was firstly used in Asia, they were called Asian options. There are two basic types: average price option and average strike option. The payoff of an average price option is the difference between the average price SAand strike price K, as follows:

Callavg.price= Max(0,(SA−K))

Put_avg.price= Max(0,(K−SA))

Average strike options, on the other hand, have a payoff based on the difference between the asset price at maturity ST and its average price SA. That is:

Call_avg.strike= Max(0,(ST−SA))

Put_avg.strike= Max(0,(SA−ST))

Asian options also differ in the type of averaging function used to calculate the payoff. We can have arithmetic Asian options or geometric Asian options.

SA(arithmetic) = 1 n " _n

∑

∏

I chose to apply the Ninomiya and Ninomiya (2008) implementation of Heston’s model on pricing Asian arithmetic average-strike options, because they are used widely for cover-ing currency exposures and protect against fluctuations of commodity prices.

Chapter 2

Theoretical framework

2.1

Notation

Symbol Explanation

t Current time

S(t) The stock price at time t

σ2(t) The annualised volatility of the underlying stock at time t

r Risk-free interest rate

α The speed of volatility mean reversion

θ The long run mean of the volatility process

β The volatility of the volatility process

ρ The correlation coefficient between S(t)and σ2(t)

B1(t), B2(t) Two independent Brownian motions

K The strike price of an Asian call option

T The expiration date of the option contract

2.2

The Heston model

In his paper Heston (1993) proposed a stochastic-volatility model which extends the tradi-tional Black–Scholes model. Heston introduced settings like the non-lognormal distribu-tion of the asset returns, the assets volatility as a stochastic process and the mean-reverting volatility.

Zhang and Jinghong Shu (2003) have shown that Heston’s model can reduce pricing errors by 25% on average compared to Black and Scholes model. Higher accuracy is exhib-ited especially in pricing middle-term and long-term options, and is explained by adding a random process for the volatility.

This paper describes the implementation of the Heston model which is a two-factor stochastic volatility(SV) model. SV models assume that the underlying asset’s price follows a stochastic process, whose volatility itself is a stochastic process, unlike the single value for volatility in the BS model. This approach overcomes the implied volatility problem from BS model and gives much better approximation for in- and out-of-the-money options.

According to Heston’s model the price of the underlying asset follows a diffusion process with drift µ and volatility σ.

dS(t) = µS(t)dt+S(t)σ(t)dB1(t), where S(0) = S0. (2.1)

The volatility, on the other hand, follows a second diffusion process. It is similar to the Cox, Ingersoll and Ross model for the evolution of interest rates. It is non-negative and

reverts to a long-term mean value θ at rate of reversion θ. The volatility of the volatility process itself is β.

dσ2(t) =α(θ−σ(t))dt+βσ(t)dB2(t), where σ2(0) =σ02. (2.2)

Let B1(t)and B2(t)be two Brownian motions with correlation ρ, while dB1(t)and dB2(t) are their respective stochastic differentials.

Cov[dB1(t), dB2(t)] =ρ dt. (2.3)

We can express the correlation relationship (2.3) into equation (2.2), and we yield the following stochastic differential equation for the volatility diffusion process:

dσ2(t) =α(θ−σ(t))dt+βσ(t)(ρ dB1(t) + q

1−ρ2dB2(t)). (2.4)

We can summarize equations (2.1) and (2.4) as Heston’s model. dS(t) = µS(t)dt+S(t)σ(t)dB1(t), where S(0) =S0,

dσ2(t) = α(θ−σ(t))dt+βσ(t)(ρ dB1(t) + q

1−ρ2dB2(t)), where σ2(0) =σ₀2.

2.3

The new Ninomiya weak approximation scheme

A new high-order scheme was developed by Kusuoka et al. (1998, 2005), for the purpose of doing very fast approximations of diffusion processes. In a paper published in 2008 Ni-nomiya S. and NiNi-nomiya M. extended the scheme and proposed a new implementation method. They used the notion of Lie algebra brackets and suggested that the stochastic approximation operator by Kusuoka can be considered a composition of the solutions of ordinary differential equations(ODEs). The set of ODEs can be approximated using the Runge–Kutta method for ODEs. This implementation is extremely fast and efficient because of the use of Runge–Kutta approximation scheme and the suggested Romberg extrapolation which boosts the performance even further.

In this paper I apply the new SDE weak-approximation scheme to pricing Asian options of European type (exercise only at maturity). The price process for the asset is modelled according to the Heston two-factor stochastic volatility model.

Without touching on details around the derivation of the scheme, we can outline it in a few basic steps. Firstly we define the pricing process:

Y1(t, x) = x1+ Z t 0 µY1(s, x)ds+ Z t 0 Y1(s, x) q Y2(s, x)dB1(s), Y2(t, x) = x2+ Z t 0 α (θ−Y2(s, x))ds+ Z t 0 β q Y(s, x)  ρ dB1(s) + q 1−ρ2dB2(s)  ,

whereΥ1describes the asset price movement, andΥ2is the diffusion process for the

volatil-ity. Here(x1, x2) ∈ (R>0)2is the initial value for price and probability respectively at time

0. (B1(t), B2(t))is a two-dimensional Brownian motion with correlation ρ,(−1<ρ<1). We want to calculate the value of an Asian average-strike option. Since it is a path-dependent option and its payoff at maturity is max(Y3(T, x)/T−K, 0), we introduce a third

process Y3(t, x):

Y3(t, x) =

Z t

The first step in the algorithm is to express the processes from Heston’s model in Itô ’s integral form: S(t) = S0+ Z t 0 rS(s)ds + Z t 0 S(s)σ(s)dB 1_(s), σ2(t) = σ02+ Z t 0 α (θ−σ(s))ds+ Z t 0 βσ(s)(ρdB 1_{(s) +}q 1−ρ2dB2(s)), I(t) = 0+ Z t 0 S(s)ds.

Let f(S(T), σ2(T), I(T)) = max{I(T)/T−K, 0} be the payoff of the Asian call option at maturity. Therefore the price of the option is D×C = D×Ef(S(T), σ2(T), I(T)) where D is the appropriate discount factor over the time to maturity. The discount factor under continuous compounding is given by:

D =exp−rT

Then we introduce the following notation in order to transform the model to Stratonovich form:

B0(s) = s,

Y(t) = (S(t), σ2(t), I(t))0, and rewrite the Itô integral form of the Heston model as

Y(t) = Y(0) + 2

∑

where ˜V0, ˜V1, ˜V2: R3 7→R3are the following functions:

˜ V0((y1, y2, y3)0) = (ry1, α(θ−y2), y1)0, ˜ V1((y1, y2, y3)0) = (y1 √ y2, ρβ √ y2, 0)0, ˜ V2((y1, y2, y3)0) = (0, β q (1−ρ2)y2, 0)0. (2.6)

Now we can step into transforming (2.5) and (2.6) in Stratonovich form.

Y(t) = Y(0) + 2

∑

where V0, V1, V2: R3 7→R3are the following functions.

V0((y1, y2, y3)0) = (y1(r−y2/2−ρβ/4), α(θ−y2) −β2/4, y1)0, V₁((y1, y2, y3)0) = (y1 √ y2, ρβ √ y2, 0)0, V2((y1, y2, y3)0) = (0, β q (1−ρ2)y2, 0)0.

To calculate this, we used formula (6.1.3) from Øksendal (2003). In our notation this formula has the form

Vj₀ =V˜j₀−1 2 2

∑

∑

After we expressed Heston’s model in Stratonovich form as it appears in (2.7) , we can introduce the new weak approximation scheme. In Theorem 1.2 from Ninomiya and Ni-nomiya (2008) we set u = 3₄. Then we have c1 = −1₂, c2 = 3₂, R11 = 3₄, R22 = 11₄, R12 = −5₄.

Let ξ1, ξ2, ξ3, and ξ4be four independent standard normal random variables. Put

S1₁ = √ 3 2 ξ1, S 1 2 = − 5√3 6 ξ1+ √ 6 3 ξ2, S2₁ = √ 3 2 ξ3, S 2 2 = − 5√3 6 ξ3+ √ 6 3 ξ4. Then E[Si_j] =0, E[Si_jSi_j00] = R_jj0δ_ii0, as desired.

The next step is to calculate the functions W1,s = ΦΨs(Z1)and W2,s = ΦΨs(Z2), s > 0

which are needed for Corollary 1.1 from Ninomiya and Ninomiya (2008). By definition of Zj, we have Z1 = − 1 2v0+ √ 3 2 ξ1v1+ √ 3 2 ξ3v2, Z2 = 3 2v0− 5√3 6 ξ1v1+ √ 6 3 ξ2v1− 5√3 6 ξ3v2+ √ 6 3 ξ4v2,

where v0, v1, and v2are the letters of an alphabet. By definition of Ψs, all terms with v0 in

the right hand side must be multiplied by s, while all the other terms must be multiplied by √ s. So, Ψs(Z1) = − 1 2sv0+ √ 3 2 ξ1 √ sv1+ √ 3 2 ξ3 √ sv2, Ψs(Z2) = 3 2sv0− 5√3 6 ξ1 √ sv1+ √ 6 3 ξ2 √ sv1−5 √ 3 6 ξ3 √ sv2+ √ 6 3 ξ4 √ sv2.

By definition ofΦ, the letters must be changed to the corresponding functions. So,

W_1,s = −1 2sV0+ √ 3 2 ξ1 √ sV1+ √ 3 2 ξ3 √ sV2, W2,s = 3 2sV0− 5√3 6 ξ1 √ sV1+ √ 6 3 ξ2 √ sV1− 5√3 6 ξ3 √ sV2+ √ 6 3 ξ4 √ sV2.

Now let’s consider the following two systems of ordinary differential equations (ODE): dy(t)

dt =Wj,s(y(t)), y(0) = y0, j=1, 2, s>0.

These systems are the ODE-valued random variables in the sense of Ninomiya and Ni-nomiya (2008), because the definition of W1,s involves random variables. Let aik be an

M×M matrix with aik = 0 for i ≤ k, and let bi be a vector of length M satisfying the

mth order conditions (46) from Ninomiya and Ninomiya (2008). The M-stage Runge–Kutta method of order m for the above mentioned systems is written as follows.

y(i) =y0+s M

∑

∑

In this case we have M = 9 and m = 7. The a and b coefficients for the seventh-order Runge-Kutta approximation scheme are taken from (Butcher, 2003) as follows:

a21 = 1 6, a32 = 1 3, a41 = 1 8, a43 = 3 8, a51 = 148 1331, a53 = 150 1331, a54 = − 56 1331, a61 = − 404 243, a63 = − 170 27 , a64 = 4024 1701, a65 = 10648 1701 , a71 = 2466 2401, a73 = 1242 343 , a74 = −19176 16807, a75 = − 51909 16807, a76 = 1053 2401, a81 = 5 154, a84 = 96 539, a85 = − 1815 20384, a86 = − 405 2464, a87 = 49 1144, a91 = − 113 32 , a93 = − 195 22 , a94 = 32 7 , a95 = 29403 3584 , a96 = −729 512, a97 = 1029 1408, a98 = 21 16, aij =0 otherwise, b =  0 0 0 32 105 1771561 6289920 243 1560 16807 74880 77 1440 11 70  .

Let n be a positive integer. Define recursively the second order scheme as the following sequence of random vectors:

Y(₀n) =Y(0),

Y(₍n_k+)_1)T/n =g(W_1,1) ◦g(W_2,1)(Y(_kT/nn) ), 0≤k ≤n−1, and define the third order scheme by the Romberg extrapolation:

X(_T2n) = 4 3Y (2n) T − 1 3Y (n) T .

Ninomiya and Ninomiya (2008) proved that the average of X(_T2n) approximates system (2.7). In particular,

lim

n→∞Ef(X

(2n)

T ) = Ef(Y(T)).

Chapter 3

Implementation of the algorithm

3.1

Monte Carlo integration

Monte Carlo integration is a numerical method for approximating the value of a definite integral. Unlike other numerical methods which use a regular grid, Monte Carlo chooses randomly the points at which the integrand is calculated.

Let’s illustrate the method with a simple example. If you want to calculate a one-dimensional integral of f(x)on the region[a, b], this means that you are interested in the area of the graph captured between a, b, f(x)and the x-axis.

Figure 3.1: Integration of f(x)over[a, b]

Using numerical methods we can approximate this area with the sum:

I = Z b a f(x)dx= (b−a) N N

∑

Under Monte Carlo integration we sample N uniformly-distributed random points on the interval [a, b] and calculate the value of the integrand at these points. Then we average the values from all sample points in order to find the area under the graph. The same ap-proach can be used for calculating multidimensional integrals. The expression for Monte Carlo approximation of the multidimensional integral over the n-dimensional unit hyper-cube is given by:

Z 1 0 · · · Z 1 0 g(x1, . . . , xn)dx1, . . . , dxn ≈ 1 N N

∑

Monte Carlo simulation has proved itself very useful for approximating the solutions of problems in the areas of mathematics, physics, biology and finance. Simulation can be considered a case of integral valuation. In order to calculate an expected value we need to integrate the probability density function of the set of outcomes. This method for simula-tion requires vectors of random numbers as an input. These vectors can be generated by a computer using deterministic methods, i.e. pseudo-random number generators.

Monte Carlo simulation has become so popular because of its characteristics:

• Simplicity of implementation and speed of calculation. In order to simulate very com-plex processes we need only to sample a set of random points x and to calculate the integrand f(x)at these points

• The inputs for the simulation (random vectors) can be generated extremely fast with modern congruential algorithms.

• The ability to approximate integrals in high dimensions without losing accuracy. The standard error does not depend on the number of dimensions, which is the main prob-lem of other methods for numerical integration.

On the other hand, Monte Carlo has certain limitations, which make it far from perfect. • This simulation method converges quite slowly. The integration error from MC is

proportional to N−12. In order to decrease the integration error 10 times, we need to

in-crease the sample size 100 times. This makes Monte Carlo simulation computationally intensive.

• The result from this method highly depends on the initial seed and the random vectors • The simulation outcome is statistical by nature, i.e. there are only probabilistic error

bounds for the estimate.

3.2

quasi-Monte Carlo integration

Quasi-Monte Carlo(QMC) integration is based on the same idea as the original Monte Carlo method. QMC, however, does not use pseudo-random numbers as an input, but sequences of quasi-random numbers with more uniform behaviour. These sequences are computer generated in a deterministic manner in order to minimize the integration error. QMC inte-gration, unlike crude Monte Carlo, is heavily dependent on the dimensions of integration. For MC the integration error is proportional to N−12, while for QMC that is (ln N)

s

N (where s

is the number of dimensions). This implies that QMC integration can achieve the same in-tegration error with considerably less sample points than crude MC. Another advantage of QMC over traditional MC integration is the nature of the error. Because quasi-Monte Carlo is based on strictly deterministic samples, we get an exact prediction for the error bound.

The quasi-random sequences used for QMC integration are also called low-discrepancy sequences. Discrepancy measures how much a given sequence deviates from the ideal se-quence which is evenly distributed in predefined geometric subset. Therefore low-discrepancy sequences fill up the unit hypercube much more evenly than a pseudo-random sequence and avoid gaps and clusterings in space.

(a) 2-dimensional pseudo-random sequence generated with java.Random

(b) the first 2 dimensions from a 4-dimensional Nieder-reiter - Xing sequence

Figure 3.2: Scatter-plot of pseudo- and quasi-random sequences, N=100

The sequence in a) does not fill up the unit space evenly, while the one in b) inhibits the grid quite uniformly.

For the pseduo-random sequence I have used Java’s built-in congruential random gen-erator and have supplied a very good seed. You can see that this sequence does not fill-ip the grid evenly. The quasi-random sequence was created with Niederreiter-Xing generator. Please notice, that the sequences on Figure 3.2 are not the same as the one used in the actual simulation. I have presented 2 uniformly-distributed sequences for better visual perception of the reader. In the simulation process I use a standard normal sequence obtained through Beasley-Springer-Moro inversion from a Niederreiter-Xing uniform sequence.

3.2.1

Low-discrepancy sequences

There are various algorithms for generating low-discrepancy sequences. Some of them are Halton, Faure, Sobol sequences. Halton is the most basic multi-dimensional quasi-random sequence. It consists of s one-dimensional sequences created with the van der Corput method using different base for every dimension. As the number of dimensions increases, it takes longer time to fill up the hypercube, which makes Halton’s quasi-random generator rather slow.

The Sobol sequence generator uses 2 as a base for all dimensions, which makes it ex-tremely fast even in high dimensions. This generator produces sequences with very good uniformity properties, and has been preferred by many researches as the best-performing algorithm. However, it requires a predefined sets of directional numbers and is not easy to implement. That is why I chose another algorithm for the purposes of my applet.

The Niederreiter-Xing random number generator was developed initially by Harold Niederreiter, and then was further developed by him and Chaoping Xing. This algorithm is based on the theory of(t, s)-sequences using the digital method. Other generators based on the digital method are Faure and Sobol. For my applet I have used Michael Meyer’s implementation of the Niederreiter-Xing algorithm. I have borrowed his java class from

the martingale package and have subsequently licensed the derivative work under the GNU General Public License. This implementation is of base 2 for all dimensions, uses the Gray code counter and bit-wise operations, which makes the algorithm extremely fast. I have also used a generator matrix by Gottlieb Pirsic for generating a 4-dimensional low-discrepancy sequence.

3.2.2

Normally distributed quasi-random sequences

The low-discrepancy generators described in the previous section produce a sequence that fills up the multidimensional unit cube evenly. Each dimension is a uniformly distributed sequence. For the purposes of financial simulations, however, we want the corresponding quasi-random numbers for the standard normal distribution. The cumulative function for the standard normal distribution is given by:

Y(x) = √1 2Π Z x −∞e −t2 2dt

Therefore we need to find the inverse cumulative distribution function and to obtain a standard normal sample from our uniform sample. The most-used method for inversion is the so-called Box-Muller algorithm. It uses a pair of uniformly distributed samples to obtain one standard normal sample. Unfortunately we cannot use it, because the Box-Muller inversion is not suitable for low-discrepancy sequences because it scrambles the sequence’s uniformity (Moro, 1995). Instead, we will use the Beasley-Springer-Moro(BSM) method. It was developed by Beasley and Springer(1977) and optimized for low-discrepancy sequences by (Moro, 1995).

Moro divided the domain of the uniform sample U in 2 regions.

• The central region of the distribution|U| ≤ 0.42 is modelled according to the Beasley-Springer algorithm: Φ−1_{(x) =U}∑3i=0aiU2i ∑4 i=0biU2i i ai bi 0 2.50662823884 1 1 -18.61500062529 -8.47351093090 2 41.39119773534 23.08336743743 3 -25.44106049637 -21.06224101826 4 3.13082909833

• For the tails of the distribution|U| > 0.42 he used truncated Chebyshev series:

Φ−1_{(x) = {} 8

∑

∑

i ai bi 0 0.3374754822726147 1 1 0.9761690190917186 0.4179886424926431 2 0.1607979714918209 4.2454686881376569 3 0.0276438810333863 4 0.0038405729373609 5 0.0003951896511919 6 0.0000321767881768 7 0.0000002888167364 8 0.0000003960315187

I have tested the inversion algorithm and evaluated the quality of the obtained stan-dard normal sequence. Then I compared this sequence with one uniform pseudo-random sequence inverted with the same method and one standard normal sequence that has not been inverted. As you can see on Figure 3.3, the inverted Niederreiter-Xing sequence is distributed in a radially-balanced way around the point (0, 0), which is the mean in 2 di-mensions.

On the table below you can also see that the sequence on the figure a) has the best charac-teristics among all three. It is therefore closest to the theoretically expected standard normal sequence.

(a) Inversion of an uniform low-discrepancy NX sequence.

(b) Inversion of an uniform pseudo-random sequence

(c) A standard normal sequence from java.Random

Figure 3.3: Inverted uniform and native standard normal sequences, N =100

Theoretical NX + Moro pseudo + Moro Java.Random(Std. Normal)

Mean 0 0.003050 0.015620 -0.017425

Median 0 0.000167 0.004800 -0.028149

Deviation 1 0.997250 1.038468 0.976299

Skewness 0 0.016785 0.073501 -0.110768

3.3

The Java applet

In this chapter I will present the Java applet for pricing Asian options. There are 3 distinct panel areas: the input panel, the control panel and the output panel. On pressing the Calcu-late button the applet starts calculating the option price, and the bar below the Reset button shows the progress of the quasi-Monte Carlo algorithm. When all simulations are complete, the option price will appear in the bottom of the Output panel, together with a graph of the price during the simulation process. The Reset button will clear the last option price from the Output panel and will reset all input variables in the Input and Controls panels to their default values. The progress bar will also be reset to 0.

Figure 3.4: The applet interface

3.3.1

Input

The input panel contains fields for entering strike price, stock price at time zero, risk-free interest rate and maturity. The last four fields in the Input panel are specific for the Heston model. According to this model the volatility of the asster price is determined by a mean-reverting stochastic process described in equation (2.2). There α is the speed of volatility mean reversion, θ is the long-run mean of the volatility and β is the volatility of the volatility

process. The last parameter in this panel ρ is the correlation of the two Brownian motions in equations (2.1) and (2.2).

In the Controls panel there are 2 additional parameters, 2 buttons and a progress bar. The number of simulations variable is crucial for the accuracy of the quasi-Monte Carlo algorithm. More simulation cycles will yield smaller integration error, and the expected option price will get closer and closer to the true value. The second variable in this panel is the number of time-steps, which sets the number of averaging days during the life of the option. For example, if time to maturity is 1 year and the user inputs 365 time-steps, then the payoff of the option will be based on the average of 365 daily prices.

Figure 3.5: Check for correctness of the correlation coefficient

(a) Message for a negative input parameter (b) Message for a non-integer input parameter

Figure 3.6: Check for correct inputs

As Ninomiya and Ninomiya (2008) points out, α, θ should be positive coefficients, and 2αθ−β2 > 0 should hold to ensure that the stochastic differential equation has an unique solution. If this inequality is not fulfilled, the user will see the following error message after pressing the Calculate button.

Figure 3.7: Check for the inequality 2αθ−β2 >0

3.3.2

Output

In the output panel of the applet there is a graph that demonstrates the convergence of the quasi-Monte Carlo method. On the Y-axis is the option price, and on the X-axis is the

num-ber of simulations performed to obtain this price. The user can change the numnum-ber of sim-ulations in the Controls panel in order to achieve the desired accuracy for the option price. The user can see graphically the speed of convergence of the quasi-Monte Carlo method in order to choose the optimal number of simlations for the desired accuracy.

Chapter 4

Analysis

4.1

Experiments

In this section I will present some experiments I did with the Java applet. I tried to inves-tigate the relationship between the price of the Asian call option and the input parameters of the Heston model. I plotted the option price over a range for each variable, keeping all other inputs constant. In this way we can see the net effect of this variable on the price.

Figure 4.1: Asian call option value depending on stock price, given K = 1.05, r = 0.05, σ =0.09, T=1, α=2, β=0.1, θ =0.09, ρ=0, time-steps = 48 and N =3000 simulations.

Figure 4.2: Asian call option value depending on strike price, given S =1, r =0.05, σ=0.09, T =1, α=2, β =0.1, θ =0.09, ρ=0, time-steps = 48 and N =3000 simulations.

Figure 4.3: Asian call option value depending on risk-free interest, given S = 1, K = 1.05, σ =0.09, T=1, α=2, β=0.1, θ =0.09, ρ=0, time-steps = 48 and N =3000 simulations.

Figure 4.4: Asian call option value depending on annualized volatility, given S = 1, K = 1.05, r = 0.05, T = 1, α = 2, β = 0.1, θ = 0.09, ρ = 0, time-steps = 48 and N = 3000 simulations.

Figure 4.5: Asian call option value depending on time to maturity, given S = 1, K = 1.05, r =0.05, σ =0.09, α=2, β=0.1, θ =0.09, ρ =0, time-steps = 48 and N =3000 simulations.

Figure 4.6: Asian call option value depending on speed of mean-reversion of the volatility, given S = 1, K = 1.05, r =0.05, σ = 0.09, T = 1, β = 0.1, θ = 0.09, ρ = 0, time-steps = 48 and N =3000 simulations.

Figure 4.7: Asian call option value depending on the long-run mean value for volatility, given S = 1, K = 1.05, r= 0.05, σ =0.09, T =1, α =2, β =0.1, ρ =0, time-steps = 48 and N =3000 simulations.

Figure 4.8: Asian call option value depending on the volatility of volatility, given S = 1, K = 1.05, r = 0.05, σ = 0.09, T = 1, α = 2, θ = 0.09, ρ = 0, time-steps = 48 and N = 3000 simulations.

Figure 4.9: Asian call option value depending on correlation between stock price and volatil-ity, given S=1, K=1.05, r =0.05, σ =0.09, T =1, α=2, β =0.1, θ =0.09, time-steps = 48 and N =3000 simulations.

Figure 4.10: Asian call option value depending on number of time-steps, given S = 1, K = 1.05, r =0.05, σ=0.09, T=1, α=2, β =0.1, θ =0.09, ρ=0 and N =3000 simulations.

4.2

Convergence and integration error

In order to evaluate the efficiency and speed of convergence of the algorithm, I measured CPU-times and compared the resulting prices. All input parameters are fixed to the default values (S=1, K =1.05, r =0.05, σ =0.09, T=1, α =2, β =0.1, θ =0.09, ρ=0, time-steps = 48), while increasing the number of simulations. The table below shows the CPU time used for calculation, the price of the option, and the absolute error. Notice that I have used the result from N=800000 simulations as a control value to measure error. The tests were performed on Intel Core2 Duo processor 2.66GHz.

Simulations CPU-time(sec) Option Price Error*

2×102 0.2 0.087305008807 0.00043424469

2×103 1.96 0.086914090857 0.00004332674

2×104 18.9 0.086874999062 0.000004234945

2×105 183.6 0.086871089883 0.000000325766

8×105 721 0.086870764117 0

Table 4.1: Performance of the applet

A quick look at the first and the last column of the table can give us an approximate measure of the speed of convergence. When the number of simulations N is increased 10 times, the absolute error decreases 10 times. Therefore we can state that the integration error of the algorithm is proportional to N−1. This is a much better measure than the crude-Monte Carlo integration error proportional to N−12 (see Section 3.1).

4.3

Accuracy of the algorithm

In the previous section I described the accuracy of the algorithm in terms of convergence and integration error. Unfortunately I did not find another Option price calculator for Asian op-tions under the Heston model, so I can not compare the correctness of the produced prices. This remains to be a subject of future investigations. It is important to point out that the Hes-ton stochastic volatility method produces more accurate results than other implied volatility methods because of its volatility mean-reverting process.

Chapter 5

Conclusion

The purpose of this thesis work, creating an applet for pricing Asian arithmetic call op-tions under the Heston stochastic volatility model using Ninomiya’s weak approximation scheme, was fully achieved. The applet implements the new approximation scheme for stochastic differential equations described in (Ninomiya and Ninomiya, 2008) and success-fully simulates the price of an Asian option. For the simulation algorithm I used the quasi-Monte Carlo method together with a low-discrepancy sequence generated by Niederreiter-Xing algorithm. The partial differential equations for the weak-approximation scheme were calculated by a 7-th order Runge-Kutta method, and Romberg extrapolation was used to boost the method to order three.

In Chapter 4 I demonstrated how the price of an Asian option changes as we change each of the input parameters for the Heston model. I also analysed the speed of calculation and the rate of convergence for the quasi-Monte Carlo method.

However there is no general conclusion for the accuracy of the pricing algorithm. This remains to be investigated further in order to compare the results from the new approxima-tion scheme with other methods incorporating the Heston diffusion processes.

Butcher, J. C. 2003. Numerical Methods for Ordinary Differential Equations, John Wiley & Sons, Chichester. Heston, S. 1993. A closed-form solution for options with stochastic volatility with applications to bond and currency options, Rev. Financ. Stud. 6, no. 2, 327-343, DOI 10.1093/rfs/6.2.327.

Karoui, N. E., M. Jeanblanc-Picque, and S. E. Shreve. 1998. Robustness of the Black and Scholes Formula, Technical Report 2.

Krykova, I. 2003. Evaluating of Path-Dependent Securities With Low Discrepancy Methods. M.S. Thesis, Worcester Polytechnic Institute.

Mikhailov, S. and U. Nögel. 2003. Heston’s Stochastic Volatilty Model: Implementation, Calibration, and some Ex-tensions., Willmot. July:74-79.

Moro, B. 1995. The Full Monte, Risk 8, february, 57-58.

Niederreiter, H. and C. Xing. 1995. Low-discrepancy sequences obtained from algebraic function fields over finite fields, Acta Arithmetica, 72, 281-298.

Ninomiya, M. and S. Ninomiya. 2008. A new weak approximation scheme of stochastic differential equations by using the Runge–Kutta method. arXiv:0709.2434v2 [math.PR] 5Jan 2008.

Øksendal, B. 2003. Stochastic Differential Equations: an Introduction with Applications, 6th ed., Springer, Berlin. Zhang, J. E. and Jinghong Shu. 2003. Pricing S&P 500 index options with Heston’s model, International Conference on Computational Intelligence for Financial Engineering, 20-23 March 2003., 85-92, DOI 10.1109/CIFER.2003.1196246, (to appear in print).

Source-code listing of Hestonwa.java

import j a v a.awt.event. ∗ ;

import j a v a.t e x t . ∗ ;

import j a v a x.swing. ∗ ;

import j a v a x.swing.border.T i t l e d B o r d e r;

import org. j f r e e .c h a r t.p l o t. ∗ ;

import org. j f r e e .c h a r t.C h a r t F a c t o r y;

import org. j f r e e .c h a r t.J F r e e C h a r t;

import org. j f r e e .data.xy.XYSeries;

11 import org. j f r e e .data.xy.X Y S e r i e s C o l l e c t i o n;

import org. j f r e e .c h a r t.ChartPanel;

/∗

∗ Java Applet f o r p r i c i n g Asian o p t i o n s

∗ under th e Heston s t o c h a s t i c v o l a t i l i t y model ∗ by Ninomiya weak approximation

∗ Version 0 . 9 8 c ∗ 2 1 . 0 4 . 2 0 0 8 ∗

∗ Author : Boyko V a s i l e v ( boyko . vasilev@gmail . com )

21 ∗

∗ WARANTY NOTICE AND COPYRIGHT

∗ This program i s f r e e s o f t w a r e ; you can r e d i s t r i b u t e i t and/or ∗ modify i t under th e terms o f th e GNU General P u b l i c L i c e n s e ∗ as published by t he Free Software Foundation ; e i t h e r v e r s i o n 2 ∗ o f t he License , or ( a t your option ) any l a t e r v e r s i o n .

∗ ∗

∗ This program i s d i s t r i b u t e d i n t he hope t h a t i t w i l l be u s e f u l , ∗ but WITHOUT ANY WARRANTY; without even the implied warranty o f

31 ∗ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE . See t h e ∗ GNU General P u b l i c L i c e n s e f o r more d e t a i l s .

∗/

p u b l i c c l a s s Hestonwa extends JApplet implements A c t i o n L i s t e n e r ,

F o c u s L i s t e n e r

{

p r i v a t e s t a t i c f i n a l long s e r i a l V e r s i o n U I D = −791580869497094349L;

// v a r i a b l e s // p a n e l s

41 p r i v a t e J P a n e l mainPanel = n u l l; p r i v a t e J P a n e l i n p u t P a n e l = n u l l; p r i v a t e J P a n e l l a b e l P a n e l = n u l l; p r i v a t e J P a n e l d a t a P a n e l = n u l l ; p r i v a t e J P a n e l c o n t r o l P a n e l = n u l l; p r i v a t e J P a n e l outputPanel = n u l l; p r i v a t e ChartPanel graphPanel1 = n u l l; // b u t t o n s p r i v a t e J B u t t o n c a l c u l a t e B u t t o n = n u l l; p r i v a t e J B u t t o n r e s e t B u t t o n = n u l l; 51 // t e x t f i e l d s p r i v a t e J T e x t F i e l d s t r i k e _ K F i e l d = n u l l; p r i v a t e J T e x t F i e l d s t o c k p r i c e _ F i e l d = n u l l ; p r i v a t e J T e x t F i e l d i n t e r e s t _ F i e l d = n u l l; p r i v a t e J T e x t F i e l d v o l a t i l i t y _ F i e l d = n u l l ; p r i v a t e J T e x t F i e l d m a t u r i t y _ F i e l d = n u l l; p r i v a t e J T e x t F i e l d a l f a _ F i e l d = n u l l; p r i v a t e J T e x t F i e l d b e t a _ F i e l d = n u l l; p r i v a t e J T e x t F i e l d t h e t a _ F i e l d = n u l l; p r i v a t e J T e x t F i e l d r h o _ F i e l d = n u l l; 61 p r i v a t e J T e x t F i e l d s i m _ F i e l d = n u l l; p r i v a t e J T e x t F i e l d dim_Field = n u l l; p r i v a t e J T e x t F i e l d p r i c e _ F i e l d = n u l l; // S t r i n g c o n s t a n t s p r i v a t e f i n a l S t r i n g CALCULATE = " C a l c u l a t e "; p r i v a t e f i n a l S t r i n g INPUT = " Input "; p r i v a t e f i n a l S t r i n g RESET = " R e s e t "; // T e x t s o f l a b e l s p r i v a t e f i n a l S t r i n g STRIKE_K_LABEL = " S t r i k e p r i c e "; p r i v a t e f i n a l S t r i n g STOCK_LABEL = " S t o c k p r i c e "; 71 p r i v a t e f i n a l S t r i n g INTEREST_LABEL = " I n t e r e s t r a t e "; p r i v a t e f i n a l S t r i n g VOLATILITY_LABEL = " V o l a t i l i t y "; p r i v a t e f i n a l S t r i n g MATURITY_LABEL = " Maturity ";

p r i v a t e f i n a l S t r i n g ALFA_LABEL = "Mean r e v e r s i o n c o e f . ( "+’ \u03B1 ’+" ) ";

p r i v a t e f i n a l S t r i n g BETA_LABEL = " V o l a t i l i t y o f Vol . ( "+’ \u03B2 ’+" ) ";

p r i v a t e f i n a l S t r i n g THETA_LABEL = " Long−term average Vol . ( "+’ \u03B8 ’ +" ) ";

p r i v a t e f i n a l S t r i n g RHO_LABEL = " C o r r e l a t i o n ( "+’ \u03C1 ’+" ) ";

p r i v a t e f i n a l S t r i n g SIM_LABEL = " S i m u l a t i o n s ";

p r i v a t e f i n a l S t r i n g DIM_LABEL = " Time−s t e p s ";

p r i v a t e f i n a l S t r i n g GRAPH1_TITLE = " Simulated o p t i o n p r i c e "; 81 p r i v a t e f i n a l S t r i n g Y1_LABEL = " Option Value ";

p r i v a t e f i n a l S t r i n g CONTROL = " C o n t r o l s ";

p r i v a t e f i n a l S t r i n g PRICE_LABEL = " Simulated o p t i o n p r i c e ";

// T o o l t i p s

p r i v a t e f i n a l S t r i n g STRIKE_K_TOOLTIP = " The p r i c e a t which one i s prepare t o e x r c i s e t h e o p t i o n "; p r i v a t e f i n a l S t r i n g STOCKPRICE_TOOLTIP = " I n i t i a l s t o c k p r i c e "; p r i v a t e f i n a l S t r i n g INTEREST_TOOLTIP = " Anual r i s k f r e e i n t e r e s t r a t e "; p r i v a t e f i n a l S t r i n g VOLATILITY_TOOLTIP = " Annualized v o l a t i l i t y o f t h e s t o c k p r i c e ";

p r i v a t e f i n a l S t r i n g MATURITY_TOOLTIP = " Maturity i n y e a r s ";

p r i v a t e f i n a l S t r i n g ALFA_TOOLTIP = "Mean r e v e r s i o n c o e f . ( "+’ \u03B1 ’+ " ) ";

91 p r i v a t e f i n a l S t r i n g BETA_TOOLTIP = " V o l a t i l i t y o f t h e v o l a t i l i t y p r o c e s s ";

p r i v a t e f i n a l S t r i n g THETA_TOOLTIP = " Long−run mean o f t h e v o l a t i l i t y p r o c e s s ";

p r i v a t e f i n a l S t r i n g RHO_TOOLTIP = " C o r r e l a t i o n ";

p r i v a t e f i n a l S t r i n g SIM_TOOLTIP = "Number o f QMC s i m u l a t i o n s ";

p r i v a t e f i n a l S t r i n g DIM_TOOLTIP = " Time−d i s c r e t i z a t i o n ";

p r i v a t e f i n a l S t r i n g PRICE_TOOLTIP = " quasi−Monte−Carlo sim ulate d p r i c e "; // E r r o r messages p r i v a t e f i n a l S t r i n g NOT_A_NUMBER = " E n t e r a number "; p r i v a t e f i n a l S t r i n g NOT_INTEGER = " E n t e r an i n t e g e r number "; p r i v a t e f i n a l S t r i n g NON_POSITIVE = " E n t e r a p o s i t i v e number "; 101 p r i v a t e f i n a l S t r i n g NON_UNIT = " E n t e r a number i n t h e i n t e r v a l [−1 ; 1 ] "; //numerical c o n s t a n t s p r i v a t e f i n a l double STRIKE_K = 1 . 0 5 ; p r i v a t e f i n a l double STOCKPRICE = 1 ; p r i v a t e f i n a l double INTEREST = 0 . 0 5 ; p r i v a t e f i n a l double VOLATILITY = 0 . 0 9 ; p r i v a t e f i n a l double MATURITY = 1 ; p r i v a t e f i n a l double ALFA = 2 ; p r i v a t e f i n a l double BETA = 0 . 1 ; p r i v a t e f i n a l double THETA = 0 . 0 9 ; 111 p r i v a t e f i n a l double RHO = 0 ; p r i v a t e f i n a l i n t SIM = 4 8 ; // number o f s i m u l a t i o n s / t r a j e c t o r i e s p r i v a t e f i n a l i n t DIM = 4 8 ; // time−d i s r e t i z a t i o n // numerical v a r i a b l e s p r i v a t e double s t r i k e _ K = STRIKE_K; p r i v a t e double s t o c k p r i c e = STOCKPRICE; p r i v a t e double i n t e r e s t = INTEREST; p r i v a t e double v o l a t i l i t y = VOLATILITY; p r i v a t e double m a t u r i t y = MATURITY; p r i v a t e double a l f a = ALFA; 121 p r i v a t e double b e t a = BETA; p r i v a t e double t h e t a = THETA;

p r i v a t e double rho = RHO;

p r i v a t e i n t sim = SIM;

p r i v a t e i n t dim = DIM;

p r i v a t e double MC;

p u b l i c i n t mc;

p r i v a t e GridBagLayout g b l = new GridBagLayout( ) ;

p r i v a t e J P r o g r e s s B a r p r o g r e s s B a r = n u l l; 131 p r i v a t e DecimalFormat numberFormatter = n u l l;

p r i v a t e double[ ] [ ] Y = new double[ 9 ] [ 3 ] ;

p r i v a t e double[ ] [ ] W = new double[ 9 ] [ 3 ] ;

p r i v a t e double[ ] yOld = new double[ 3 ] ;

NX nxs = n u l l;

Derivs f S i n g l e = n u l l;

Derivs fDouble = n u l l;

141 p u b l i c void i n i t ( ) {

// g e t c o n t e n t pane

Container contentPane = getContentPane( ) ;

// c r e a t e main panel

mainPanel = new J P a n e l( ) ;

// s e t box l a y o u t

mainPanel.s e t L a y o u t(new BoxLayout(mainPanel,BoxLayout.Y_AXIS) ) ;

// add main panel t o c o n t e n t pane

contentPane.add(mainPanel) ;

151

// S e t proper number format

DecimalFormatSymbols symbols = new DecimalFormatSymbols( ) ;

symbols.s e t D e c i m a l S e p a r a t o r( ’ . ’) ;

numberFormatter = new DecimalFormat(" # # # . # # # # # # # # # # # # ",symbols) ;

// c r e a t e in pu t panel

i n p u t P a n e l = new J P a n e l(new GridLayout( 0 , 3 ) ) ;

i n p u t P a n e l.s e t B o r d e r(new T i t l e d B o r d e r(INPUT) ) ; 161 // add i t t o t h e main panel

mainPanel.add(i n p u t P a n e l) ;

// c r e a t e l a b e l P a n e l

l a b e l P a n e l = new J P a n e l(new GridLayout( 0 , 1 ) ) ;

i n p u t P a n e l.add(l a b e l P a n e l) ;

// c r e a t e data panel

d a t a P a n e l = new J P a n e l(new GridLayout( 0 , 1 ) ) ;

// add i t t o in pu t panel

171 i n p u t P a n e l.add(d a t a P a n e l) ;

// c r e a t e C o n t r o l panel i n s i d e In pu t Pa ne l

c o n t r o l P a n e l = new J P a n e l(new GridLayout( 7 , 2 ) ) ;

c o n t r o l P a n e l.s e t B o r d e r(new T i t l e d B o r d e r(CONTROL) ) ;

i n p u t P a n e l.add(c o n t r o l P a n e l) ;

// c r e a t e output panel

outputPanel = new J P a n e l(g b l) ;

outputPanel.s e t B o r d e r(new T i t l e d B o r d e r(" Output ") ) ; 181

// add i t t o t h e main panel

mainPanel.add(outputPanel ,BorderLayout.CENTER) ;

// add l a b e l s

J L a b e l l a b e l = new J L a b e l(STRIKE_K_LABEL) ;

l a b e l P a n e l.add(l a b e l) ;

s t r i k e _ K F i e l d = new J T e x t F i e l d( ) ; 191 // add t o o l t i p

s t r i k e _ K F i e l d .s e t T o o l T i p T e x t(STRIKE_K_TOOLTIP) ;

// s e t value

s t r i k e _ K F i e l d .s e t T e x t(numberFormatter.format(STRIKE_K) ) ;

// add f o c u s l i s t e n e r

s t r i k e _ K F i e l d .a d d F o c u s L i s t e n e r(t h i s) ;

// add i t t o data panel

d a t a P a n e l.add(s t r i k e _ K F i e l d) ; // c r e a t e s t o c k l a b e l 201 l a b e l = new J L a b e l(STOCK_LABEL) ; l a b e l P a n e l.add(l a b e l) ; // c r e a t e s t o c k f i e l d s t o c k p r i c e _ F i e l d = new J T e x t F i e l d( ) ; // add t o o l t i p s t o c k p r i c e _ F i e l d.s e t T o o l T i p T e x t(STOCKPRICE_TOOLTIP) ; // s e t value

s t o c k p r i c e _ F i e l d.s e t T e x t(numberFormatter.format(STOCKPRICE) ) ;

// add f o c u s l i s t e n e r

s t o c k p r i c e _ F i e l d.a d d F o c u s L i s t e n e r( t h i s) ; 211 // add i t t o data panel

d a t a P a n e l.add(s t o c k p r i c e _ F i e l d) ; // c r e a t e i n t e r e s t l a b e l l a b e l = new J L a b e l(INTEREST_LABEL) ; l a b e l P a n e l.add(l a b e l) ; // c r e a t e s t o c k f i e l d i n t e r e s t _ F i e l d = new J T e x t F i e l d( ) ; // add t o o l t i p 221 i n t e r e s t _ F i e l d .s e t T o o l T i p T e x t(INTEREST_TOOLTIP) ;

// s e t v a l u e P a n e l . add ( dataPanel , BorderLayout . CENTER) ;

i n t e r e s t _ F i e l d .s e t T e x t(numberFormatter.format(INTEREST) ) ;

// add f o c u s l i s t e n e r

i n t e r e s t _ F i e l d .a d d F o c u s L i s t e n e r(t h i s) ;

// add i t t o data panel

d a t a P a n e l.add(i n t e r e s t _ F i e l d) ; // c r e a t e v o l a t i l i t y l a b e l 231 l a b e l = new J L a b e l(VOLATILITY_LABEL) ; l a b e l P a n e l.add(l a b e l) ; // c r e a t e v o l a t i l i t y f i e l d v o l a t i l i t y _ F i e l d = new J T e x t F i e l d( ) ; // add t o o l t i p v o l a t i l i t y _ F i e l d.s e t T o o l T i p T e x t(VOLATILITY_TOOLTIP) ; // s e t value

v o l a t i l i t y _ F i e l d.s e t T e x t(numberFormatter.format(VOLATILITY) ) ;

// add f o c u s l i s t e n e r

// add i t t o data panel d a t a P a n e l.add( v o l a t i l i t y _ F i e l d) ; // c r e a t e m a t u r i t y l a b e l l a b e l = new J L a b e l(MATURITY_LABEL) ; l a b e l P a n e l.add(l a b e l) ; // c r e a t e m a t u r i t y f i e l d 251 m a t u r i t y _ F i e l d = new J T e x t F i e l d( ) ; // add t o o l t i p m a t u r i t y _ F i e l d.s e t T o o l T i p T e x t(MATURITY_TOOLTIP) ; // s e t value

m a t u r i t y _ F i e l d.s e t T e x t(numberFormatter.format(MATURITY) ) ;

// add f o c u s l i s t e n e r

m a t u r i t y _ F i e l d.a d d F o c u s L i s t e n e r(t h i s) ;

// add i t t o data panel

d a t a P a n e l.add(m a t u r i t y _ F i e l d) ; 261 // c r e a t e a l f a l a b e l l a b e l = new J L a b e l(ALFA_LABEL) ; l a b e l P a n e l.add(l a b e l) ; // c r e a t e a l f a f i e l d a l f a _ F i e l d = new J T e x t F i e l d( ) ; // add t o o l t i p a l f a _ F i e l d.s e t T o o l T i p T e x t(ALFA_TOOLTIP) ; // s e t value

271 a l f a _ F i e l d.s e t T e x t(numberFormatter.format(ALFA) ) ;

// add f o c u s l i s t e n e r

a l f a _ F i e l d.a d d F o c u s L i s t e n e r(t h i s) ;

// add i t t o data panel

d a t a P a n e l.add(a l f a _ F i e l d) ; // c r e a t e b e t a l a b e l l a b e l = new J L a b e l(BETA_LABEL) ; l a b e l P a n e l.add(l a b e l) ; 281 // c r e a t e b e t a f i e l d b e t a _ F i e l d = new J T e x t F i e l d( ) ; // add t o o l t i p b e t a _ F i e l d.s e t T o o l T i p T e x t(BETA_TOOLTIP) ; // s e t value

b e t a _ F i e l d.s e t T e x t(numberFormatter.format(BETA) ) ;

// add f o c u s l i s t e n e r

b e t a _ F i e l d.a d d F o c u s L i s t e n e r(t h i s) ; 291 // add i t t o data panel

d a t a P a n e l.add(b e t a _ F i e l d) ;

// c r e a t e t h e t a l a b e l

l a b e l P a n e l.add(l a b e l) ; // c r e a t e t h e t a f i e l d t h e t a _ F i e l d = new J T e x t F i e l d( ) ; // add t o o l t i p 301 t h e t a _ F i e l d.s e t T o o l T i p T e x t(THETA_TOOLTIP) ; // s e t value

t h e t a _ F i e l d.s e t T e x t(numberFormatter.format(THETA) ) ;

// add f o c u s l i s t e n e r

t h e t a _ F i e l d.a d d F o c u s L i s t e n e r( t h i s) ;

// add i t t o data panel

d a t a P a n e l.add(t h e t a _ F i e l d) ; // c r e a t e rho l a b e l 311 l a b e l = new J L a b e l(RHO_LABEL) ; l a b e l P a n e l.add(l a b e l) ; // c r e a t e rho f i e l d r h o _ F i e l d = new J T e x t F i e l d( ) ; // add t o o l t i p r h o _ F i e l d.s e t T o o l T i p T e x t(RHO_TOOLTIP) ; // s e t value

r h o _ F i e l d. s e t T e x t(numberFormatter.format(RHO) ) ;

// add f o c u s l i s t e n e r

321 r h o _ F i e l d.a d d F o c u s L i s t e n e r(t h i s) ;

// add i t t o data panel

d a t a P a n e l.add(r h o _ F i e l d) ; // c r e a t e sim l a b e l l a b e l = new J L a b e l(SIM_LABEL) ; c o n t r o l P a n e l.add(l a b e l) ; // c r e a t e sim f i e l d 331 s i m _ F i e l d = new J T e x t F i e l d( ) ; // add t o o l t i p s i m _ F i e l d.s e t T o o l T i p T e x t(SIM_TOOLTIP) ; // s e t value

s i m _ F i e l d.s e t T e x t(numberFormatter.format(SIM) ) ;

// add f o c u s l i s t e n e r

s i m _ F i e l d.a d d F o c u s L i s t e n e r(t h i s) ;

// add i t t o data panel

c o n t r o l P a n e l.add(s i m _ F i e l d) ; 341 // c r e a t e dim l a b e l l a b e l = new J L a b e l(DIM_LABEL) ; c o n t r o l P a n e l.add(l a b e l) ; // c r e a t e dim f i e l d dim_Field = new J T e x t F i e l d( ) ; // add t o o l t i p

dim_Field.s e t T o o l T i p T e x t(DIM_TOOLTIP) ;

// s e t value

351 dim_Field. s e t T e x t(numberFormatter.format(DIM) ) ;

// add f o c u s l i s t e n e r

dim_Field.a d d F o c u s L i s t e n e r(t h i s) ;

// add i t t o data panel

c o n t r o l P a n e l.add(dim_Field) ;

// c r e a t e c a l c u l a t e button

c a l c u l a t e B u t t o n = new J B u t t o n(CALCULATE) ;

// add a c t i o n l i s t e n e r

361 c a l c u l a t e B u t t o n.a d d A c t i o n L i s t e n e r( t h i s) ;

// add i t t o button panel

c o n t r o l P a n e l.add(c a l c u l a t e B u t t o n) ;

// c r e a t e r e s e t button

r e s e t B u t t o n = new J B u t t o n(RESET) ;

// add a c t i o n l i s t e n e r

r e s e t B u t t o n.a d d A c t i o n L i s t e n e r(t h i s) ;

// add i t t o button panel

c o n t r o l P a n e l.add(r e s e t B u t t o n) ; 371 G r i d B a g C o n s t r a i n t s c = new G r i d B a g C o n s t r a i n t s( ) ; // c r e a t e p r o g r e s s bar p r o g r e s s B a r = new J P r o g r e s s B a r( ) ; p r o g r e s s B a r.s e t S t r i n g P a i n t e d(t r u e) ; p r o g r e s s B a r.s e t V a l u e( 0 ) ; p r o g r e s s B a r.s e t I n d e t e r m i n a t e( f a l s e) ; //add i t t o p r o g r e s s panel c. f i l l = G r i d B a g C o n s t r a i n t s.HORIZONTAL; c.gridwidth = 3 ; 381 c o n t r o l P a n e l.add(progressBar ,c) ; J F r e e C h a r t c h a r t; // G r i d B a g C o n s t r a i n t s c = new G r i d B a g C o n s t r a i n t s ( ) ; c. f i l l = G r i d B a g C o n s t r a i n t s.BOTH; // c r e a t e panel f o r J f r e e c h a r t

graphPanel1 = new ChartPanel(n u l l) ;

c.ipady = 0 ; // r e s e t t o d e f a u l t

391 c.weightx = 1 . 0 ; // r e q u e s t any e x t r a h o r i z o n t a l space

c.weighty = 0 . 9 2 ; // r e q u e s t any e x t r a v e r t i c a l space

c.anchor = G r i d B a g C o n s t r a i n t s.PAGE_END; //bottom o f space

c.i n s e t s = new I n s e t s( 0 , 0 , 0 , 0 ) ; //top padding

c.g r i d x = 0 ;

c.gridy = 0 ;

c.gridwidth = 2 ;

outputPanel.add(graphPanel1,c) ;

// c r e a t e p r i c e l a b e l

c.gridwidth = 1 ; c.g r i d x = 0 ; c.gridy = 1 ; c.weighty = 0 . 0 8 ; l a b e l 2.s e t S i z e( 0 , 2 0 ) ; outputPanel.add(l a b e l 2 ,c) ; // c r e a t e p r i c e f i e l d p r i c e _ F i e l d = new J T e x t F i e l d( ) ; 411 // add t o o l t i p p r i c e _ F i e l d.s e t T o o l T i p T e x t(PRICE_TOOLTIP) ; // s e t value p r i c e _ F i e l d.s e t T e x t(numberFormatter.format( 0 ) ) ; p r i c e _ F i e l d.s e t E d i t a b l e( f a l s e) ;

// add i t t o data panel

c.g r i d x = 1 ; c.gridy = 1 ; c.weighty = 0 . 0 8 ; outputPanel.add(p r i c e _ F i e l d ,c) ; 421 outputPanel.v a l i d a t e ( ) ; } // end i n i t

p u b l i c void focusGained(FocusEvent e) { }

p u b l i c void f o c u s L o s t(FocusEvent e) { // who c a l l e d f o c u s L o s t ? O b j e c t s o u r c e = e.g e t S o u r c e( ) ; // i f s t r i k e K f i e l d 431 i f (s o u r c e == s t r i k e _ K F i e l d) { // check s t r i k e _ K s t r i k e _ K = r e a d P o s i t i v e(s t r i k e _ K F i e l d , s t r i k e _ K , s t r i k e _ K F i e l d.g e t T o o l T i p T e x t( ) ) ; r e t u r n; } // I f s t o c k p r i c e f i e l d i f (s o u r c e == s t o c k p r i c e _ F i e l d) { // check s t o c k p r i c e 441 s t o c k p r i c e = r e a d P o s i t i v e(s t o c k p r i c e _ F i e l d , s t o c k p r i c e , s t o c k p r i c e _ F i e l d.g e t T o o l T i p T e x t( ) ) ; r e t u r n; } // I f i n t e r e s t r a t e f i e l d i f (s o u r c e == i n t e r e s t _ F i e l d) { // check i n t e r e s t i n t e r e s t = r e a d P o s i t i v e(i n t e r e s t _ F i e l d , i n t e r e s t , 451 i n t e r e s t _ F i e l d .g e t T o o l T i p T e x t( ) ) ; r e t u r n; }

// I f v o l a t i l i t y f i e l d i f (s o u r c e == v o l a t i l i t y _ F i e l d) { // check v o l a t i l i t y v o l a t i l i t y = r e a d P o s i t i v e(v o l a t i l i t y _ F i e l d , v o l a t i l i t y , v o l a t i l i t y _ F i e l d .g e t T o o l T i p T e x t( ) ) ; 461 r e t u r n; } // I f m a t u r i t y f i e l d i f (s o u r c e == m a t u r i t y _ F i e l d) { // check m a t u r i t y m a t u r i t y = r e a d P o s i t i v e(m a t u r i t y _ F i e l d , maturity , m a t u r i t y _ F i e l d.g e t T o o l T i p T e x t( ) ) ; r e t u r n; 471 } // I f a l f a f i e l d i f (s o u r c e == a l f a _ F i e l d) { a l f a = r e a d P o s i t i v e(a l f a _ F i e l d , a l f a , a l f a _ F i e l d .g e t T o o l T i p T e x t( ) ) ; r e t u r n; } 481 // I f b e t a f i e l d i f (s o u r c e == b e t a _ F i e l d) { b e t a = readDouble(b e t a _ F i e l d , beta , b e t a _ F i e l d.g e t T o o l T i p T e x t( ) ) ; r e t u r n; } // I f t h e t a f i e l d i f (s o u r c e == t h e t a _ F i e l d) { 491 t h e t a = r e a d P o s i t i v e(t h e t a _ F i e l d , t h e t a , t h e t a _ F i e l d .g e t T o o l T i p T e x t( ) ) ; r e t u r n; } // I f rho f i e l d i f (s o u r c e == r h o _ F i e l d) { rho = readUnit(r h o _ F i e l d , rho, 501 r h o _ F i e l d.g e t T o o l T i p T e x t( ) ) ; r e t u r n; } // I f SIM f i e l d i f (s o u r c e == s i m _ F i e l d) { sim = r e a d I n t(sim_Field ,

sim, s i m _ F i e l d.g e t T o o l T i p T e x t( ) ) ; r e t u r n; 511 } // I f DIM f i e l d i f (s o u r c e == dim_Field) { dim = r e a d I n t(dim_Field , dim, dim_Field.g e t T o o l T i p T e x t( ) ) ; r e t u r n; } } 521

p u b l i c void actionPe rformed(ActionEvent e) {

// determine , who c a l l e d a c t i o n l i s t e n e r

O b j e c t s o u r c e = e.g e t S o u r c e( ) ;

// I f C a l c u l a t e button

i f (s o u r c e == c a l c u l a t e B u t t o n) {

// v e r i f y c o n d i t i o n f o r unique s o l u t i o n o f SDEs and c o r r e l a t i o n

i f ( ( 2 ∗a l f a∗t h e t a−b e t a∗beta>0) && (Math.abs(rho) <=1) ) { 531 X Y S e r i e s C o l l e c t i o n d a t a s e t = new X Y S e r i e s C o l l e c t i o n( ) ;

XYSeries p r i c e S e r i e s = new XYSeries(" P r i c e ") ;

d a t a s e t.r e m o v e A l l S e r i e s( ) ; p r i c e S e r i e s .c l e a r ( ) ; // Prepare p r o g r e s s bar p r o g r e s s B a r.setMaximum( 1 0 0 ) ; p r o g r e s s B a r.s e t V a l u e( 0 ) ; 541 long b e f o r e = System.c u r r e n t T i m e M i l l i s( ) ; double t ; i n t pb; MC= 0 . 0 ;

f S i n g l e = new Derivs(m a t u r i t y/dim, i n t e r e s t , rho, beta, a l f a , t h e t a) ;

fDouble = new Derivs(m a t u r i t y/(2∗dim) , i n t e r e s t , rho, beta , a l f a ,

t h e t a) ; nxs = new NX( 4 ) ; // i n i t N i e d e r r e i t e r−Xing g e n e r a t o r i n 4 dimensions f o r (mc= 0 ; mc<=sim; mc++) { t=GetX2n( ) ; 551 MC=MC+t; nxs. r e s t a r t( ) ;

pb=mc∗100 / sim; // r e f r e s h p r o g r e s s b a r and add p o i n t s t o c h a r t // only 100 ti me s f o r a l l s i m u l a t i o n s

i f (p r o g r e s s B a r.getValue( ) ! = pb) {

p r o g r e s s B a r.s e t V a l u e(pb) ;

i f (mc ! = 0 ) p r i c e S e r i e s.add(mc,MC/mc) ; } MC=MC/sim; 561 p r i c e _ F i e l d.s e t T e x t(numberFormatter.format(MC) ) ; J F r e e C h a r t c h a r t ; d a t a s e t.a d d S e r i e s(p r i c e S e r i e s) ; c h a r t = C h a r t F a c t o r y.createXYLineChart(GRAPH1_TITLE, " I t e r a t i o n s " , Y1_LABEL, d a t a s e t , P l o t O r i e n t a t i o n.VERTICAL, 571 tr ue , f a l s e , f a l s e) ; graphPanel1.s e t C h a r t(c h a r t) ; } e l s e {

i f ( ( 2 ∗a l f a∗t h e t a−b e t a∗beta<=0) ) JOptionPane.showMessageDialog

(n u l l ,

" P l e a s e observe 2 "+’ \u03B1 ’+’ \u03C1 ’+"−"+’ \u03B2 ’+" ^2>0 i n order t o ensure unique s o l u t i o n t o Heston ’ s

SDEs " ,

" I n v a l i d i n p u t s " ,

581 JOptionPane.ERROR_MESSAGE) ;

i f (Math.abs(rho) >1) JOptionPane.showMessageDialog(n u l l , " C o r r e l a t i o n "+’ \u03C1 ’+" should be i n t h e i n t e r v a l [−1 ; 1 ] " , " I n v a l i d i n p u t s " , JOptionPane.ERROR_MESSAGE) ; } } 591 // I f R e s e t button i f (s o u r c e == r e s e t B u t t o n) { // r e s e t a l l T e x t F i l e d s and v a r i a b l e s t o t h e i n i t i a l v a l u e s s t r i k e _ K = STRIKE_K;

s t r i k e _ K F i e l d.s e t T e x t(numberFormatter.format(STRIKE_K) ) ;

s t o c k p r i c e = STOCKPRICE;

s t o c k p r i c e _ F i e l d.s e t T e x t(numberFormatter.format(STOCKPRICE) ) ;

i n t e r e s t = INTEREST;

i n t e r e s t _ F i e l d .s e t T e x t(numberFormatter.format(INTEREST) ) ;

601 v o l a t i l i t y = VOLATILITY;

v o l a t i l i t y _ F i e l d.s e t T e x t(numberFormatter.format(VOLATILITY) ) ;

m a t u r i t y = MATURITY;

m a t u r i t y _ F i e l d.s e t T e x t(numberFormatter.format(MATURITY) ) ;

a l f a = ALFA;

b e t a = BETA;

b e t a _ F i e l d.s e t T e x t(numberFormatter.format(BETA) ) ;

t h e t a = THETA;

t h e t a _ F i e l d .s e t T e x t(numberFormatter.format(THETA) ) ;

611 rho = RHO;

r h o _ F i e l d.s e t T e x t(numberFormatter.format(RHO) ) ;

sim = SIM;

s i m _ F i e l d.s e t T e x t(numberFormatter.format(SIM) ) ;

dim = DIM;

dim_Field.s e t T e x t(numberFormatter.format(DIM) ) ;

p r i c e _ F i e l d .s e t T e x t(numberFormatter.format( 0 ) ) ; p r o g r e s s B a r.s e t V a l u e( 0 ) ; r e t u r n; } 621 } p r i v a t e double readUnit(J T e x t F i e l d f i e l d , double oldValue, S t r i n g t i t l e ) { boolean isOK = t r u e; double newValue = 1 ; t r y { // t e s t in pu t

newValue = Double.parseDouble(f i e l d .g e t T e x t( ) ) ;

631 }

c a t c h (NumberFormatException e) { // ERROR message

JOptionPane.showMessageDialog(n u l l , NOT_A_NUMBER, t i t l e , JOptionPane.ERROR_MESSAGE) ; isOK = f a l s e ; }

i f (Math.abs(newValue) >1) { // ERROR message

JOptionPane.showMessageDialog(n u l l , 641 NON_UNIT, t i t l e , JOptionPane.ERROR_MESSAGE) ; isOK = f a l s e; } i f (isOK) { r e t u r n newValue; } e l s e {

f i e l d.s e t T e x t(numberFormatter.format(oldValue) ) ; 651 r e t u r n oldValue; } } p r i v a t e double readDouble(J T e x t F i e l d f i e l d , double oldValue, S t r i n g t i t l e ) { boolean isOK = t r u e; double newValue = 1 ;

t r y { // t e s t in pu t

661 newValue = Double.parseDouble(f i e l d.g e t T e x t( ) ) ;

}

c a t c h (NumberFormatException e) { // ERROR message

JOptionPane.showMessageDialog(n u l l , NOT_A_NUMBER, t i t l e , JOptionPane.ERROR_MESSAGE) ; isOK = f a l s e; } i f (isOK) { 671 r e t u r n newValue; } e l s e {

f i e l d. s e t T e x t(numberFormatter.format(oldValue) ) ;

r e t u r n oldValue; } } p r i v a t e double r e a d P o s i t i v e(J T e x t F i e l d f i e l d , double oldValue, 681 S t r i n g t i t l e ) { boolean isOK = t r u e; double newValue = 1 ; t r y { // t e s t in pu t

newValue = Double.parseDouble( f i e l d.g e t T e x t( ) ) ; }

c a t c h (NumberFormatException e) { // ERROR message

JOptionPane.showMessageDialog(n u l l , NOT_A_NUMBER, t i t l e , 691 JOptionPane.ERROR_MESSAGE) ; isOK = f a l s e; }

i f (newValue <= 0 ) { // ERROR message

JOptionPane.showMessageDialog(n u l l , NON_POSITIVE, t i t l e , JOptionPane.ERROR_MESSAGE) ; isOK = f a l s e; } 701 i f (isOK) { r e t u r n newValue; } e l s e {

f i e l d. s e t T e x t(numberFormatter.format(oldValue) ) ;

r e t u r n oldValue; } } // Reads i n t e g e r numbers 711 p r i v a t e i n t r e a d I n t(J T e x t F i e l d f i e l d , i n t oldValue,

S t r i n g t i t l e) { boolean isOK = t r u e; i n t newValue = 1 ; t r y { // t e s t in pu t newValue = I n t e g e r.p a r s e I n t(f i e l d.g e t T e x t( ) ) ; }

c a t c h (NumberFormatException e) { // ERROR message

JOptionPane.showMessageDialog(n u l l , 721 NOT_INTEGER, t i t l e , JOptionPane.ERROR_MESSAGE) ; isOK = f a l s e; }

i f (newValue <= 0 ) { // ERROR message

JOptionPane.showMessageDialog(n u l l , NON_POSITIVE, t i t l e , JOptionPane.ERROR_MESSAGE) ; 731 isOK = f a l s e; } i f (isOK) { r e t u r n newValue; } e l s e {

f i e l d. s e t T e x t(numberFormatter.format(oldValue) ) ;

r e t u r n oldValue; }

} 741

double moro_inv(double u) {

/∗ r e t u r n s t h e i n v e r s e o f cumulative normal d i s t r i b u t i o n f u n c t i o n Reference > The F u l l Monte , by B o r i s Moro , RISK 1 9 9 5 ( 2 ) ∗/

double a[ ] = { 2 . 5 0 6 6 2 8 2 3 8 8 4 , −18.61500062529 , 4 1 . 3 9 1 1 9 7 7 3 5 3 4 , −2 5 . 4 4 1 0 6 0 4 9 6 3 7 } ; double b[ ] = {−8 . 4 7 3 5 1 0 9 3 0 9 0 , 751 2 3 . 0 8 3 3 6 7 4 3 7 4 3 , −21.06224101826 , 3 . 1 3 0 8 2 9 0 9 8 3 3 } ; double c[ ] = { 0 . 3 3 7 4 7 5 4 8 2 2 7 2 6 1 4 7 , 0 . 9 7 6 1 6 9 0 1 9 0 9 1 7 1 8 6 , 0 . 1 6 0 7 9 7 9 7 1 4 9 1 8 2 0 9 , 0 . 0 2 7 6 4 3 8 8 1 0 3 3 3 8 6 3 , 0 . 0 0 3 8 4 0 5 7 2 9 3 7 3 6 0 9 , 0 . 0 0 0 3 9 5 1 8 9 6 5 1 1 9 1 9 , 0 . 0 0 0 0 3 2 1 7 6 7 8 8 1 7 6 8 , 761 0 . 0 0 0 0 0 0 2 8 8 8 1 6 7 3 6 4 , 0 . 0 0 0 0 0 0 3 9 6 0 3 1 5 1 8 7 } ; double x,r; x=u−0 . 5 ; i f (Math.abs(x) < 0 . 4 2 ) {