extfieldsofstudy:
Economics
Computer Science
Mathematics

text:
Branching Particle Pricers with Heston Examples

The use of sequential Monte Carlo within simulation for path-dependent option pricing is proposed and evaluated. Recently, it was shown that explicit solutions and importance sampling are valuable for efficient simulation of spot price and volatility, especially for purposes of path-dependent option pricing. The resulting simulation algorithm is an analog to the weighted particle filtering algorithm that might be improved by resampling or branching. Indeed, some branching algorithms are shown herein to improve pricing performance substantially while some resampling algorithms are shown to be less suitable in certain cases. A historical property is given and explained as the distinguishing feature between the sequential Monte Carlo algorithms that work on path-dependent option pricing and those that do not. In particular, it is recommended to use the so-called effective particle branching algorithm within importance-sampling Monte Carlo methods for path-dependent option pricing. All recommendations are based upon numeric comparison of option pricing problems in the Heston model.

Introduction
Several years ago, Heston (1993) introduced an asset model that has withstood the test of time and become widely popular. Its popularity stems from the availability of a closed-form solutions for European option prices and the inclusion of stochastic volatility. It is, in a way, the "Black-Scholes model of the stochastic volatility world" and is widely used to model stock, bond and foreign currency prices. Recently, Kouritzin (2018) introduced explicit weak solutions to the two-dimensional stochastic differential equation (SDE) solutions of the Heston model that are easy to simulate.
Option prices that do not allow closed-form solutions are often priced by Monte Carlo simulation. Using Kouritzin (2018)'s explicit solutions, it is possible to simulate multiple Heston paths {(S j t , V j t ), t ≥ 0} N j=1 efficiently. For reasons that will become apparent in the next paragraphs, we often refer to each of these simulated paths as particles. Each particle (S j t , V j t ) has the desired distribution (the one under which we want to price) under the pricing measure P , and the associated likelihood (or Radon-Nikodym derivative) L j t = P P j Ft can be seen as its weight. For this reason, the simulation algorithm for option pricing of Kouritzin (2018) can be thought of as the analog of a weighted particle filter in sequential Monte Carlo methods, producing a weighted empirical measure of the form σ N [0,t] at each time step. Here, δ (S j [0,t] ,V j [0,t] ) denotes Dirac measure and (S j [0,t] , V j [0,t] ) denotes the path of the j th -particle's price and volatility over [0, t] but held constant from t on. The purpose of the present paper is to discuss the use of resampling and branching in the option pricing procedure in much the same manner as resampling and branching are used in sequential Monte Carlo methods.
The weighted particle filter (credited to Handschin (1970), Handschin & Mayne (1969)) also produces an object of the form of σ N [0,t] as an approximation of the unnormalized filter in problems like tracking and model selection. While the application and the exact likelihood form differ, one can think of the overall object in the same way in the context of Monte Carlo Heston simulation. The particle likelihoods {L j t } N j=1 provide relative odds that the simulated paths {(S j [0,t] , V j [0,t] )} N j=1 are good representation of the paths (S [0,t] , V [0,t] ) of the underlying model. However, there is a well known problem in particle filtering (or sequential Monte Carlo): If left alone, as t increases, many of the simulated particles {(S i t , V i t )} N i=1 move away from the underlying model (S t , V t ), effectively shrinking the particle system to fewer and fewer relevant particles (i.e. particles that are good representations of the model).
The solution in sequential Monte Carlo was to resample or branch the particles regularly in an unbiased manner so as to kill bad particles while replicating the good ones with high probability. The first and still most popular algorithm, called the bootstrap algorithm, introduced by Gordon et al. (1993), simply redistributes the different particles independently back to their last observed locations according to their relative likelihoods. This algorithm was later improved by the residual resampling of Liu & Chen (1998), stratified resampling of Kitagawa (1996), combined residual-stratified resampling discussed in Douc & Cappé (2005) and systematic resampling of Carpenter et al. (1999). In all cases, enough resampling was done that the weights L i t could be reset to 1 after resampling. A different approach was taken in Kouritzin (2017b). Here, particles are branched individually rather than resampled collectively and branching only occurs when there is sufficient need. Indeed, the cost of rebalancing the system by resampling and branching is the introduction of noise that immediately degrades estimates. Too much resampling or branching is thus undesirable. With this limited branching, weights are kept after branching. Different branching methods, the residual, combined, dynamic and effective particle branching algorithms, were presented and shown to outperform all of the resampled particle algorithms mentioned above on tracking and model selection (see Kouritzin (2017b)). The simplest of these algorithms, the residual algorithm, has also been analyzed in Kouritzin (2017a).
Herein, we focus on two of these new branching algorithms and introduce them in the Heston simulation framework of Kouritzin (2018). Separately, we show empirically that the bootstrap algorithm is less suitable for pricing options that are strongly path-dependent, mainly because its excessive resampling affects the distribution of path estimates. Conversely, we demonstrate easy deployment and effectiveness of the new branching algorithms and recommend the use of the so-called effective particle branching within simulation-based, path-dependent option pricing problems when ν = nκ 2 4 for any n ∈ N (for there is no need to consider branching or resampling when ν = nκ 2 4 ). To do so, we focus on pricing path-dependent options and options allowing for early exercise.
Path-dependent options, including American, Asian and swap options, can be effectively priced using simulation and dynamic programming (DP). Paths of the underlying option prices and other relevant market variables (e.g. stochastic volatility) are first simulated up to maturity of the option. Then, the DP algorithms run backwards in time through each simulation (i.e. particle) comparing the payoff resulting from immediate exercise to the discounted future (risk-neutral) expected payout if the option is not exercised. In particular, American (and other continuouslyexecutable) options are approximated discretely, resulting in Bermudan-style options whose optimal exercise time is determined recursively starting at the option maturity date. A key breakthrough in this approach by Longstaff & Schwartz (2001) and others (see Carriere (1996), Tsitsiklis and Van Roy (2001)) was the realization that this conditional expectation of discounted future payout could be estimated using the simulated particle paths and, further, that the (uncomputable) conditional expectations could be approximated by easily-computable projections. Initially, least-squares regression was used in this projection process. However, this requires solving an often-ill-conditioned system of linear equations. Kouritzin (2018) then proposed replacing regression with one-step stochastic approximation (SA), which can allow more efficient, accurate pricing without numerical stability issues, but whose performance is sensitive to the chosen step gain. To reduce the sensitivity of the algorithm to the selection of the step gain, we recommend a modification of the algorithm based on Polyak and Juditsky (1992). Herein, we further develop this Monte Carlo pricing approach, started by Longstaff & Schwartz (2001), by introducing importance sampling into the approach. Finally, we demonstrate the advantage of resampling through the use of branching algorithms over resampling with interacting algorithms such as the multinomial bootstrap.
The rest of this paper is laid out as follows: In the next section, we recall the Heston model's weak solutions and present the weighted Heston simulation algorithm in a framework suggestive of resampling or branching. In Section 3, we recall the bootstrap algorithm and explain why it is not suitable for American-style option pricing. We also introduce the branching algorithms and insert them into the weighted Heston algorithm to yield new combined and effective particle branching Heston simulators. In Section 4, we insert these branching simulators into the SA/DP option pricing algorithm. Section 5 consists of our empirical results regarding option pricing. Finally, Section 6 contains our conclusions.

Explicit and Weighted Heston Algorithms
Heston (1993) developed a closed form solution for European call options in a specific stochastic volatility model. Broadie & Kaya (2006) then came up with an exact simulation method for the Heston model. Both of these methods are based upon characterizing and evolving the distribution of (S t , V t ). Later, Kouritzin (2018) showed that the Heston SDEs were actually explicitly solvable in a way that is very convenient for simulation, which facilitates the pricing of exotic options whose prices cannot be written in closed-form. Specifically, if we let for some n ∈ N, then the following result is developed in Kouritzin (2018): Theorem 1. Let ε ∈ (0, 1); T > 0; V 0 , S 0 be given random variables with V 0 > ε; and {W 1 , . . . , W n , B} be independent standard Brownian motions on (Ω, F , {F t } t∈[0,T ] , P ). Let also Then, η ε is a stopping time and L t∧ηε is a likelihood (i.e. a non-negative L r -martingale with respect to P for any r > 0 such that E[ L t∧ηε ] = 1 for all t ≥ 0). Moreover, (B, β) are independent standard Brownian motions and If we take n = 4ν κ 2 + 1 2 ∨ 1 in (2.1) (or equivalently in (2.3)), then is called the closest explicit Heston model. Then, Theorem 1 provides a means to produce weighted particles of the desired Heston model until the volatility drops below ε and to resort to the closest explicit model thereafter.
Remark 1 (Notation). We are using S, V for solutions to the closest explicit Heston model under P , reserving S, V for the solution to (1.1). Henceforth, we will use It is necessary to stop (at η ε ) before the volatility gets too small. Indeed, when the volatility reaches V 1 2 t = 0, the closest explicit and general Heston volatility equations become deterministic d V t = ν κ dt, dV t = νdt and it is obvious which solution one has. The model distributions are then singular to each other when ν k = ν.
copies of ( L, S, V , η ε ), where S, V , L, η ε are as in Theorem 1, and S j [0,T ] , V j [0,T ] denotes the (continuous) paths of the j th particle. Then, by the (measure-valued) strong law of large numbers, we find that where the right hand side means the distribution of (S [0,T ] , V [0,T ] ) from (2.5). However, the weights { L j T ∧ηε } N j=1 will typically be very uneven so the effective number of particles, approximated here by the ratio N ef f , is much smaller than N and the convergence rate is slow (in N). We will introduce resampling and branching to prevent this from happening.

Sequential Monte Carlo and Branching Heston Algorithm
The results of the previous sections can be applied to the efficient simulation of Heston paths and their associated likelihood, which are then used to weight each path in option price calculations. In order to avoid the case in which a small number of paths have significantly higher likelihoods than the rest, thus increasing the variance of the estimate and reducing the performance of the simulation algorithm, we resort to regular resampling of the particles.
We first recall the bootstrap algorithm, which will later be shown to be less suitable for American option pricing. We then present two branching algorithms that can improve the performance of the Heston simulation procedure of Kouritzin (2018).
3.1. Bootstrap Algorithm. The main idea behind the bootstrap algorithm is to construct a multinomial distribution using the likelihoods associated with each particle and to sample from this distribution. For clarity, in Algorithm 1, we summarize the bootstrap algorithm as it might be used within the simulation of the Heston model stock prices and volatilities.
Remark 2. The step described in line 5 was purposely vague so as not to detract from the bootstrap algorithm. These Heston simulation details appear below and can be inserted in place of this step. However, we argue and show experimentally that this bootstrap algorithm does not perform as well as other resampling algorithms for American-type option pricing and we recommend not using it.
The problem with the bootstrap algorithm in path-dependent option pricing is that there is too much resampling noise. Suppose that we kept track of each particle through all times. When resampling occurs, a "child particle" moves to a "parent" particle site and one can create historical paths {S j [0,T ] , V j [0,T ] } N j=1 by appending the child's path to its parent's. If we just focus on discrete times, then it is shown in Del Moral et al. (2001) that for the bootstrap algorithm, meaning the limit is not the (discrete-time) process distribution of the desired Heston model but rather product measure of the marginals, i.e. independence. This means that the historical property (2.6) does not hold and estimating the conditional expectation of future discounted payoff can not be determined correctly from the cross-sectional data when the bootstrap is used, which has a significant impact on path-dependent option pricing. If one mistakenly calculated this estimate with bootstrap, he or she would wind up with a regular expectation instead of a conditional one, which is not useful to many option pricing problems. A good way to correct for this is to switch to the branching algorithms in Kouritzin (2017b). 3: for t = 1 to T do Evolve the particles independently 5:

Normalize weights
7: for j = 1 to N do 8: end for 10: end for 20: end procedure We note that there is a second (non-bootstrap) resampling algorithm in Del Moral et al. (2001) that does have the required historical property. However, this second algorithm is slower and arguably more complicated than the branching algorithms discussed below.
3.2. Branching Algorithms. Compared to collective resampling algorithms such as the bootstrap one presented above, individual branching improves speed by simplifying decisions. The decision to branch (or kill) a particle (that is, the decision to resample the same particle more than once, or not at all) is based upon the weight of that individual particle and the total weight of all particles, and not on the individual values of the other ones. Our basic branching framework is described in Algorithm 2.
The key steps of Algorithm 2 are described between lines 7 and 22. They serve determine the new number of particles N t and weights {L j t } Nt j=1 in an unbiased manner. The main idea of the branching algorithm is the following. When the prior Algorithm 2 Branching Algorithm 1: procedure GeneralBranching(S 0 ,V 0 ,N,T ,r) 2: for t = 1 to T do Evolve particles independently using Theorem 1 3: Check which particles to branch 5: end for Branching part of the algorithm 14: Let p be a random permutation of {l + 1, l + 2, ..., N t−1 } 17:

20:
for k = 1 to N j t do 21: end for 25: end for 26: end procedure weight L j t for particle j is extreme (that is, outside of a certain interval around the average weight), we do (limited) residual-style branching, which helps to preserve the process distribution. Particles that are branched result in zero or more particles, which are assigned the average weight A t , are added at the same location as the parent. In other words, we copy (or kill) the path with extreme prior weight and give the copies, if there are any, the current average weight. When its prior weight L j t is not extreme, a particle is not branched and gets to keep its prior weight. The flexibility in this class of algorithms is how we determine "extreme".
The resampling parameter r = (r t ; t ≥ 0) determines the size of the interval around the average weight A t outside of which particles are considered extreme and are branched. This parameter can be obtained using different methods (see below) and controls the amount of branching: If r t = ∞, then there is no branching. If r t = 1, then every particle is branched. Both of these cases are rarely good. Instead, one looks for a good fixed r ∈ (1, ∞) (combined branching) or even makes r t depend upon the system of particles (effective particle branching).
When it is necessary to record the whole path and not just its final value, we just append a child's path to that of its parent. It follows that each child at time t ≥ 2 has a unique parent at time t − 1 and grandparent t − 2 but a particle at time t might have no or several children at time t + 1. Then, the historical paths (or a discretization thereof) can be used for option pricing. Indeed, in Section 5, we price options using The precise mathematical proof will resemble that of the strong law of large numbers in Theorem 5.1 in Kouritzin (2017a) but is left to future work. More information on this branching algorithm framework can be found in Kouritzin (2017b).
One key aspect of these branching algorithms discussed in Kouritzin (2017b) is that the number of particles do not vary wildly. Wild particle number variation could affect performance and cause a method to fail. Notice that the average weight A t in line 4 is normalized by the initial number of particles (rather than the current), which forces the expected number of future particles given the current to be the initial N.
We present below two branching options. The second one is a refinement of the first one, and it should perform better. We mention that the performance of the branching algorithms is tied to how well they control their particle numbers, as well as under what circumstances they will branch more, and direct the reader to Kouritzin (2017b) for more details and other branching algorithms.
3.2.1. Combined Branching. This is the basic algorithm used in Algorithm 2. It uses stratified resampling to help control the number of particles and improve algorithm performance. The {ρ j t } N t−1 j=1 are negatively correlated, which is desirable for particle control. The generation of a large number of new particles will be more likely to be compensated after by the generation of a smaller number of new particles, because of the negative correlation between {ρ j t } N t−1 j=1 . In lines 7 to 12 of Algorithm 2, we handle the non-branched particles while moving the ones that will be branched to the front of the line. In line 15, we create the required number of stratified uniform random numbers {U j t } and in lines 18 to 22, we branch the particles designated for branching. We have ρ j t = 1 so we are actually using a residual technique (see Kouritzin (2017b)). However, generating negatively correlated {ρ j t } N t−1 j=1 , rather than independent one as in the Residual branching algorithm of Kouritzin (2017b) reduces the probability of getting mostly large or mostly small uniform random numbers. This reduces the variation in ρ j t and thus in the number of particles. Remark 3. In combined branching, we use a suitable fixed r t = r > 1. The better values are model dependent and our choices are given within each example in Section 5. For most problems, the optimal proportion of paths that are branched at each time step is between 0.05 and 0.65. For a given problem, r can be chosen by minimizing the root-mean-square error (RMSE), or another measure of the difference, between the theoretical and estimated values of quantities that are known theoretically, such as the expectation of L T or S T . This minimization can be performed using a stochastic gradient descent algorithm. In many cases, the RMSE will not be very sensitive to r. When this is the case, one should choose r to manage the trade off between adding branching noise now and re-distributing the particles for better future use.

Effective Particle Branching.
In the prior method, a constant r t was used for simplicity only. If one is in need of fast, accurate option pricing, then there are better choices.
The uneven weights that arise from direct use of Theorem 1 result in an effective number of particles, N ef f , that can be estimated. In our setting, the estimated number of effective and non-effective particles are: From the second equality, we see that N ef f t−1 = N t−1 if all L j t are the same so all particles are equally effective and N ef f t−1 = 1 if all but one of the L j t were 0 so there is only one effective particle. (None of the L j t can be zero but they can be arbitrarily close.) Otherwise, it gives us a number somewhere in between that can be interpreted as the effective number of particles. It is very reasonable to anticipate better results when branching either more or fewer particles in the situation there are few effective ones. A first intuition might lead us to the conclusion that it is better to branch more in order to obtain more effective particles immediately. However, if those few particles with high weights happen to be wrong, then we will likely be creating a large number of copies of these "bad" paths. We do not assume either a priori but rather, in the effective particle branching algorithm, we set for experimentally determined constants c ef f , c nonef f > 0 and let the data decide.
3.3. Weighted and Explicit Heston Simulation. We now turn our attention to line 3 in the basic branching algorithm (this step also appears in line 5 of the bootstrap algorithm), in which ( j=1 . Here, we recall how to perform this step using Theorem 1 (see also Kouritzin (2018)).
Defining constants we find that (2.2,2.4) can be rewritten as which simplifies the simulation of the values at a future time step. Indeed, the stochastic integral in (3.3) is conditionally (given V ) Gaussian since V and B are independent so it can be simulated as a centered normal random variable with variance a 2 t t−1 V s ds. Even the likelihood (3.4) avoids stochastic integrals. There are a number of ways to compute the two deterministic integrals in (3.3) and (3.4). Kouritzin (2018) mentions that Simpson's 1 3 rule with M = 6 is a good choice. In general, the choice of M should depend on the Heston's model parameters, more specifically on the ratio ν κ 2 , and on the size of the time steps. Indeed, if ν κ 2 is low (or in our case, if 4ν κ 2 is close to 2) a larger M, i.e. a finer discretization, should be used in order for the simulated V j t 's not to drop too close to zero. When ν κ 2 is higher, numerical tests have shown that the precision of the estimate is not affected by the use of a smaller M (see Section 5.1). When this is the case, it is better to use a lower value for M to speed up the algorithm. Of course, if the time steps used in the simulation are longer, a larger M should be selected to increase the precision of the deterministic integral. However, since one of the advantages of Algorithm 3 is its speed, we recommend using it to price path-dependent options when the financial market needs to be simulated at short time intervals.
For the presentation of the algorithm, it will be notationally convenient to remove the tildes and to define two more constants

4̺
, and α M = e − ̺ 2M . (3.5) The resulting procedure to simulate paths of the Heston model using T time steps of length 1 is given in Algorithm 3.
Remark 4. Algorithm 3 is presented in its most basic form. Its performance can be further improved through the use of antithetic variates for the Z j,i 's. Another way to further speed up the algorithm could be to generate a large sample of Normal random variables and to re-use them in such a way to minimize the dependence thus created. The application of this method is left for future work. 3: for t = 1 to T do 4: end for 10: 12:  We now turn our attention to pricing American-style options in the Heston model. The risk-neutral (pricing) measure we use going forward is the one denoted by P in (2.5). In other words, ν = nκ 2 4 for any n ∈ N, and we must resort to Algorithm 3, along with branching, to simulate paths of (S, V ) for pricing purposes.
We consider an option that can be exercised only once, at any time before or at maturity. The payoff process of the option represents the amount received by the investor given that the option is exercised at t, for t ∈ {0, 1, . . . , T } 1 We denote this process by {Z t , t ≥ 0} and assume that it satisfies either for some non-negative measurable function p, where R t = 1 t t k=1 S t is the running average of the underlying asset price. For example, using the notation x + = max(x, 0), we have that p(t, S t ) = e −µt (K − S t ) + for an American put or p(t, R t ) = e −µt (K − R t ) + for an Asian put with early exercise. Recall from (2.5) that µ denotes the drift of S t under the risk-neutral measure, and so e −µt discounts the payoff from time t to 0.
Remark 5. In numerical implementations, the running average price R t should be recursively calculated from the price S by setting R 0 = 0 and The price at t = 0 of an American-style option with payoff process Z t , denoted here by C 0 is defined as where T t,T denotes the collection of stopping times taking values in {t, t + 1, ..., T } and E denotes expectation with respect to P . By letting τ t ∈ T t,T be such that American option pricing problems therefore involve calculating E[Z τt ] for different values t ∈ {0, . . . , T }, which can be done using Monte Carlo simulations and dynamic programming. This so-called least-squares Monte Carlo method has the advantage of being highly flexible (see e.g. Longstaff & Schwartz (2001), Clément et al. (2002)). Here, we recall the main ideas of the method and present an algorithm that makes use of stochastic approximation techniques to increase execution speed and numerical stability (see also Kouritzin (2018)).
As explained below, a key step of the pricing algorithm relies on the projection of the future payoff. Clément et al. (2002) used the general assumptions discussed below to ensure that this projection can be done correctly. We first note that the procedure relies on the underlying model being a Markov chain {X t } ∞ t=0 and on the payoff being of the form Z t = f (t, X t ). Therefore, when pricing payoffs depending on the running average, it is necessary to add a third state variable, R t , to our underlying model. We denote by D = D S × D V × D R the state space of the Markov process (S, V, R). Then, the following assumptions are required to use projections to approximate E[Z τt ].
• Total: there are measurable R-valued functions (e k ) ∞ k=1 on D such that where e J = (e 1 , ..., e J ) ′ . The point of these conditions is that we need a suitable collection of functions {e k } ∞ k=1 that can be used for projection of a future payoff. Typically, is an ordering of the basis functions Theoretically, suitable stopping times can be obtained by working backwards from T : and ∪{Z t = 0} do not affect the recursion in (4.2).
In practice, the conditional expectations { E[Z τ t+1 |F t ]} T −1 t=0 and stopping times are not immediately computable. Therefore, we follow Longstaff & Schwartz (2001) and approximate the conditional expectations by projections onto the closed linear span of e J (S t , V t ) (see also Clément et al. (2002), Kouritzin (2018)), which leads to the following approximate stopping times, which depends on the number J of projection functions used The problem is still not solved because this projection, is computed using coefficients α J t that are not known. However, Longstaff & Schwartz (2001) showed that these coefficients α J t can be estimated using cross-sections of a Monte Carlo particle system via linear regression. Kouritzin (2018) then showed that stochastic approximation can be a favorable alternative to the original least-squares approach. Neither of these works contemplated using the method in a branching particle system, which is what we do here.
In addition, to improve the robustness and the performance of the algorithm, we use the average of the coefficients α J t in the projection. Averaging the coefficients was first proposed by Polyak and Juditsky (1992) and is now a widely used method that can speed up the convergence of the stochastic approximation algorithm. We recall our branching particle system and apply the stochastic approximation algorithm to our setting next. Hereafter, suppose {(L j t , S j t , V j t ), t ≥ 0} N j=1 is as in Section 3.2 so (3.1) holds, where (S, V ) satisfies (2.5) and Z j t = p(t, S j t ) or Z j t = p(t, R j t ), with R j t = 1 t t k=1 S j t . Moreover, define the particle-by-particle stopping times and suppose is defined by α J,0 t = 0 and k = 1, and by , k + 1), Z j t > 0, for j = 1, 2, ..., N. Then, as explained in Kouritzin (2018) s. for any γ > 0 and which with help of Clément et al. (2002) yields and gives a means to estimate the option value E[Z τ 0 ] as well as optimal execution τ 0 . This procedure is described in Algorithm 4. Algorithm 4 can be improved by using the average of the coefficientsᾱ J,N . This modification is done by initializingᾱ J t = 0 for all t ∈ {0, . . . , T − 1}, by adding the following line after Line 11: /k, and by calculating the projection in line 15 using e J (S j t , V j t ) ′ᾱJ t . Generally, the exponent χ in the gain step γ k χ is equal to 1 when there is no averaging. However, the best value for χ is generally below 1 with averaging. With the right parameters, averaging the coefficients speeds up convergence of the algorithm. Another important advantage of this modification is that it significantly increases the robustness of the algorithm with respect to the parameters γ and χ; the new algorithm is thus easier to use in practice because convergence is achieved for wider intervals of parameters.
There is still is a potential problem since the projection coefficients depend upon the stopping times and the stopping times depend upon the projection coefficients. Fortunately, the dependence can be decoupled: α J,j In this section, we illustrate how branching can be used to increase the accuracy and efficiency of Algorithm 3. This algorithm was shown in Kouritzin (2018) to be an interesting tool to price path-dependent options in the Heston model. However, a disadvantage of the simulation method is that it relies on likelihood, or weights, whose variance increase in time, possibly reducing the accuracy of Monte Carlo price estimators.
Here we consider three different parameter sets to identify cases where branching clearly improves the performance of the weighted simulation algorithm. The parameters considered are given in Table 1. 7 We show that, for a given number of simulations, branching is often very effective in increasing the precision of price estimates and in decreasing their variance.
To perform the simulations, we discretize the time horizon using m time steps per year. We first consider both Asian and European straddle options with maturity T , 6 R j only need to be simulated when pricing Asian options with early exercise.

7
The parameters were chosen to illustrate the efficiency of the branching algorithms rather than for their fit to a specific set of financial data.

Algorithm 4 Stochastic Approximation
1: procedure AmericanOptionPricing(e k , γ, χ, N, T ) 2: , copies of (L,Ŝ,V , R) 6 5: if Z j t > 0 then 10:   In all cases, the use of the explicit weak solution of Theorem 1 necessitates the simulation of weights L j since 4ν κ 2 is not an integer.
We first price options without early exercise to highlight the effect of branching on the efficiency of the simulation algorithm. We then price American-type options using the algorithm of Section 4 and show that when re-sampling is needed, branching algorithms should be used instead of the bootstrap algorithm. 5.1. Efficiency of the Simulation Algorithms. In this section, we price financial options using the branching algorithms introduced in Section 3.2 in combination with the weighted simulation method (Algorithm 3). The Monte Carlo price estimates are calculated using the weighted simulated price paths and we assume that no early exercise is allowed. In other words, the price estimate of an option with payoff Z T is given by (5.1) Simulations are obtained using Algorithm 3 with Simpson's 1 3 rule with 50 time steps per year and M sub-intervals (see Table 2 for specific values) to compute the deterministic integrals. The branching parameters r t (for combined branching) and c ef f , c nef f (for effective branching), and the stopping parameter ε were chosen in the following way: (1) Find ε * minimizing 1 N N j=1 (L j T − 1) . (2) Using ε * , find r t (for combined branching) or c ef f , c nef f (for effective branch- Using the known quantities E[L T ] = 1 and E[S T ] = e µT S 0 allows for efficient calibration of the branching parameters. We note however that, depending on the problem, the objective function is not smooth, since it is based on simulations, and often almost flat around its minimum, which makes it more difficult to identify optimal parameters. However, it is possible to identify intervals of good parameters that can be used in the algorithms. We note that in the cases of PS2 and PS3, 4ν κ 2 is close to 2, which is the threshold for the Feller condition. Thus, while the Feller condition ensuring that V t does not hit 0 is satisfied in all cases, we expect V t to drop much closer to 0 when using PS2 and PS3, compared to PS1. Since the explicit solution underlying the weighted simulation algorithm only holds until V t gets arbitrarily close to 0, parameter sets that keep V t away from 0 should be "easier" to simulate from. This explains why for the second and the third parameter set, we need smaller time sub-intervals to simulate V t (that is, M = 6) and a higher ε. Indeed, setting ε = 10 −5 helps catch the few paths where the volatility drops too much, causing the weights to explode. When this happens, the weight is set to 0; the particle is effectively dropped from the Monte Carlo estimate, inducing a bias, which we show to be small. There exists situations where branching is crucial to obtaining a fast and accurate pricing algorithm. We illustrate this situation by pricing a European straddle option (for which a closed form is available) using the market parameters of Table 1 and the simulation parameters of Table 2.
That is, we calculate the Monte Carlo estimator of the option price C using (5.1).
To test the accuracy of the various pricing algorithms, we use the relative RMSE where C is the real option price calculated using the integral form available for European option prices in the Heston model (see Heston (1993) and Albrecher et al. (2007)). We compute the relative RMSE for N ∈ {10 4 , 5 × 10 4 , 10 5 , 5 × 10 5 } simulations. Figure 1 shows that, for the second and third parameter sets (PS2 and PS3), the relative RMSE of the price estimate is significantly reduced by branching. Branching is however less necessary for the first parameter set (PS1), for which the weighted algorithm 3 alone results in very precise estimates.
As mentioned above, the branching parameters r t , c ef f and c nef f were chosen to optimize the performance of the pricing algorithm. Depending on the parameters used, the resulting optimal proportion of branched particles vary. This proportion is much lower for PS1 (under 20%) than for PS2 and PS3 (between 30% and 35%).
The choice of M, the number of intervals used to approximate the numerical integrals, depends on the parameters. Table 3 shows the sensitivity of the RMSE of the price estimate to changes in M. To isolate the impact of M on the RMSE, we perform all calculations without resampling the particles. The results obtained with PS1 are already very precise with M = 2; increasing it does not reduce the RMSE, but slows down the algorithm. This is not true for PS2 and PS3; increasing M has a significant impact of the RMSE. We note that in the case of PS2, using M = 8 could lead to more precise results, while also further increasing the execution time. We also observe that for a fixed M, the algorithm is much slower for PS1, since it requires the simulation of n = 9 Gaussian processes (compared to n = 3 for PS2 and PS3). However, using M = 2 allows the run time to remain low. 5.1.2. Performance of the Branching Algorithms. We have shown that for a given number of initial particles, branching can increase the accuracy of the price estimate. This additional step can however slow down the pricing algorithm. In this section, we consider the standard deviation of the price estimate as a function of time. In addition to European options, we also price Asian options with fixed strike K = 100, for which there does not exist a closed-form expression. Because of its speed, the weighted simulation algorithm is particularly well-suited for path-dependent options (see also Kouritzin (2018)). To compare the performance of the algorithms, we consider the standard deviation of the estimate relative to the true option price, which we approximate by computing the price estimate 50 times. The true option price is computed using the semi-analytical formula in the case of the European option. The mean of 50 price estimates obtained using N = 5 × 10 5 and the weighted algorithm with effective particle branching is considered to be the "true" price of the Asian option.
There are cases in which branching clearly increases the performance of the pricing algorithm. This is true for PS2 and PS3; it can be seen in Figures 2 and 3 that the standard deviation of the price estimator can be reduced much faster when branching is used. That is, although branching does slow down the pricing procedure, the resulting improvement is such that fewer initial particles are needed for the standard deviation of the price estimate to reach a level similar to the one reached without branching for a larger number of initial particles. For example, in the case of PS2 (see Figure 2(B)), the European price estimate using either branching algorithms and N = 10 5 has a smaller standard deviation than the non-branched estimate obtained with five times more simulated paths. In this case, pricing with branching is more than 3 times faster. Significant improvements due to branching are also observed for the third parameter set and when pricing Asian options (see Figure 3).
In the case of PS1, the algorithm is already very precise without branching; additional resampling does not seem to be necessary. However, it does not significantly increase the RMSE or the standard deviation of the European option price estimate. Branching does however increase the standard deviation of the Asian option price estimate. This is most likely due to the path-dependent nature of the Asian option payoff; branching can lead to path degeneracy (the surviving particles having few distinct ancestors), which impacts the variance of the price estimate. Although path degeneracy does not appear be a significant problem in our numerical examples in general, the very high performance of the weighted simulation algorithm (without resampling) for PS1 leads to the variance of the price estimate increasing when resampling is added.
The need for resampling is linked to the size of the ratio 4ν κ 2 . Indeed, this ratio is much higher for PS1 ( 4ν κ 2 = 8.50); we could say that the Feller condition is "completely satisfied" and V t stays well away from 0. Further numerical tests also show that the variance of the likelihood L T is much lower than in the case of the other two parameter sets. PS2 and PS3 have much smaller ratios 4ν κ 2 ; they are equal to 2.65 and 2.50, respectively. For such parameters, V t has a higher probability of dropping closer to 0, which increases the variance of the likelihoods and creates the need for more branching.
Our numerical results show that in general, when Algorithm 3 is used, a branching algorithm should be implemented alongside. Effective branching, of which combined branching is a special case, should be prioritized. When necessary, depending on the market parameters considered, branching can be turned off by setting c ef f and c nef f high enough such that no particle is branched. Improvements due to branching are most significant when the variance of the weights L is high; for example, when 4ν κ 2 is close to 2 (as with PS2 and PS3) or when the option to price has a longer maturity, since the variance of the particle system tends to increase with time.
Remark 7. In our setting, the performance of the combined branching algorithm is similar to the one of the effective particle branching one. This is likely due to the fact that, in the cases we study, the number of effective particles remains relatively constant as t increases. Thus, using both algorithms results in resampling a similar proportion of the particles at each time step, whether r t is kept constant (combined branching) or not (effective particle branching). In situations where the distribution of the weights can vary widely from one time step to the other (in particle filtering, for example), the increased flexibility provided by the effective particle branching algorithm can result in a significant performance improvement (see Kouritzin (2017b)). Since it is similar to the combined branching algorithm in terms of complexity of implementation and computation ressources, we still recommend the use of the effective particle branching algorithm, even in cases where its advantage is less obvious.

5.2.
Comparison with Bootstrap Algorithm. In this section, we highlight an important advantage of the branching algorithms of Section 3.2 compared to the bootstrap algorithm of Section 3.1: the preservation of the process distribution. To do so, we price Asian options with early exercise, whose price strongly depends on the distribution of the historical paths. We focus on the second and the third sets of parameters (PS2 and PS3) since they are the ones for which the performance of the pricing algorithm is the most improved by resampling.
To price the options, we use the SA/DP algorithm exposed in Section 4 using the first 3 √ J weighted Laguerre polynomials as basis functions for each of volatility, price and average price. We also use the averaged coefficientsᾱ J,N t in order to increase the robustness of the algorithm to the choice of parameters γ and χ. In our setting, compared to the well-known least-squares Monte Carlo algorithm, SA/DP allows us to use a large number of basis functions without sacrificing numerical stability and at a lower computing cost than the least-square Monte Carlo algorithm. This leads to more precise price estimates (see Kouritzin (2018)). As explained in Section 3, the amount of noise introduced by bootstrap resampling can affect the historical law of the particle system, which is problematic when pricing options with early exercise via recursive conditional expectations determined from cross-sectional data, as in Algorithm 4. We price at-the-money Asian call options with fixed strike and early exercise (that is, we consider Bermudan-style options that can be exercised at every discrete time step). To do so, we use Algorithm 4 with Z j t = max(R j t − K, 0) and K = S 0 . We consider options with one year to maturity and use the same market and simulation parameters as in Section 5.1. Figures 4 and 5(A) present the prices obtained using two resampling methods: effective branching and bootstrap resampling (see Algorithm 1). In each case, the empirical distribution presented is estimated using 50 price estimates. These results are compared to a "ground truth price" of 7.67 for PS2 and 6.89 for PS3, indicated by the red dashed line in the figures. These prices were computed without resampling, but with a large number of simulations (N = 10 6 ) and basis functions (J = 12 3 ), to obtain a very precise result. Table 4 in the Appendix contains the SA/DP parameters used for all calculations. These parameters control the speed at which the gain step γ k χ goes to 0 as k increases in the stochastic approximation algorithm. In our case, we do not want the gain steps to decrease too fast (that is, γ should not be too small and χ should not be too large), as we want each of the N paths to have some impact on the estimates α J,N t . In addition, since we use averaged coefficientsᾱ J,N t to estimate the conditional expectations, we choose χ < 1, as recommended in Section 4. Nonetheless, if γ is too large, or if χ is too small, then we may not average out the noise fast enough and can run out of particles before getting to a good estimate forᾱ J,N t . In Kouritzin (2018), a stochastic search method was used to find reasonable γ and χ. However, the averaging method that we employ herein is far less sensitive to particular choice of γ and χ. Within the interval of values for which the gain step decreased at a satisfying rate, we ran the algorithm for different values of γ and χ and chose those that minimized 1 N N j=1 (Z j τ J,j − e J (S j t , V j t ) ′ᾱ J,N t ) 2 for t ∈ 0, . . . , T − 1. Figures 4 and 5 clearly show that branching algorithms should be preferred over multinomial bootstrap when resampling is necessary to price American-type options. Indeed, the prices obtained with effective particle branching are in general closer to the true price than when bootstrap is used. For the second parameter set (PS2), prices obtained with resampling also show a clear convergence towards the true price as J increases, while such a convergence is not observed when using bootstrap resampling. The prices obtained with the third parameter set (PS3) are not as clear in terms of convergence; a larger number of simulations is needed to reduce the bias of the price estimate.
It has been shown theoretically in Del Moral et al. (2001) that bootstrap resampling kills the process distribution of the simulated paths as the number of simulations gets very large. Our results support this empirically; the price estimates obtained using bootstrap do not converge as clearly to the true price, and the variance of the estimator is in all cases significantly higher than with effective particle branching. For this reason, branching algorithms, and effective branching in particular, should be preferred when importance sampling is used for pricing options with early exercise. Resampling algorithms based on particle branching are also typically faster than bootstrap algorithms, since less resampling occurs.

Conclusions and Future Work
Herein, we show that the simulation method presented by Kouritzin (2018) can be improved by resampling, especially when 4ν κ 2 is close to 2. Since bootstrap resampling affects the process distribution of the particle system, the branching algorithms of Kouritzin (2017b), also presented in Section 3.2 in the context of option pricing, are a better alternative. In particular, effective particle branching offers more flexibility and allows the level of resampling to be a function of the state of the system. We also highlight that the use of stochastic approximation algorithms within the dynamic programming framework for American option pricing is still valid when historical price paths are obtained via branching resampling algorithms. Our use of the average coefficientsᾱ J,N in the pricing algorithm increases the robustness of the algorithm to the choice of parameters γ and χ. This is an important improvement over the algorithm used by Kouritzin (2018), as it increases the stability of the pricing method and makes it easier to use, since less work is needed to find the optimal parameters. Future work will explore more advanced ways for the algorithms to "learn" the optimal resampling parameters. Additional improvements linked to recycling simulated random variables in order to speed up the simulation algorithm will also be considered.
Appendix A. Additional Information on Numerical Results Table 4 contains the stochastic approximation parameters used to obtain the numerical results of Section 5.2.


extfieldsofstudy:
Computer Science
Engineering
Mathematics

text:
Acoustic Scene Classification Using Bilinear Pooling on Time-liked and Frequency-liked Convolution Neural Network

The current methodology in tackling Acoustic Scene Classification (ASC) task can be described in two steps, preprocessing of the audio waveform into log-mel spectrogram and then using it as the input representation for Convolutional Neural Network (CNN). This paradigm shift occurs after DCASE 2016 where this framework model achieves the state-of-the-art result in ASC tasks on the (ESC-50) dataset and achieved an accuracy of 64.5%, which constitute to 20.5% improvement over the baseline model, and DCASE 2016 dataset with an accuracy of 90.0% (development) and 86.2% (evaluation), which constitute a 6.4% and 9% improvements with respect to the baseline system. In this paper, we explored the use of harmonic and percussive source separation (HPSS) to split the audio into harmonic audio and percussive audio, which has received popularity in the field of music information retrieval (MIR). Although works have been done in using HPSS as input representation for CNN model in ASC task, this paper further investigate the possibility on leveraging the separated harmonic component and percussive component by curating 2 CNNs which tries to understand harmonic audio and percussive audio in their natural form, one specialized in extracting deep features in time biased domain and another specialized in extracting deep features in frequency biased domain, respectively. The deep features extracted from these 2 CNNs will then be combined using bilinear pooling. Hence, presenting a two-stream time and frequency CNN architecture approach in classifying acoustic scene. The model is being evaluated on DCASE 2019 sub task 1a dataset and scored an average of 65% on development dataset, Kaggle Leadership Private and Public board.

I. INTRODUCTION
Over the past few years, there has been an increase in interest in the research on acoustic/audio analysis [8]. The possible reasons for this phenomenon can be linked to the break-through from visual classification [1] and more publicly available platform and communities for researchers or practitioner to build and benchmark their classification models and discuss their findings and architecture. Detection and Classification of Acoustic Scenes and Events (DCASE) [9] is one such platform, and acoustic scene classification has been one of the on-going challenges that appear in DCASE challenges throughout the years.
Acoustic Scene Classification (ASC) is in general, the ability to computationally differentiate various scene (e.g. beach, airport, streets, shopping centre and etc.) from a given audio signal and the potential of such research is vast (e.g. audio surveillance systems [10], robots or autonomous mobile [11] and environmental noise monitoring [12], [13]).
However, ASC is a challenging task, an audio recording from natural environment can contain multiple sound sources, ambient sounds and each sound sources can also have differing characteristics [10]. This sound sources are usually interleaving and overlapping in a single audio recording, proofing it difficult in identifying the number of sound sources and even if we can separate the sound source, recognizing the quasi-stationary and non-stationary sounds is also challenging due to their changing nature.
A popular paradigm in tackling ASC tasks can be broken down into two steps. First step which is pre-processing of the audio signal and second steps, using the feature extracted from preprocessing fit it into a classifier. The sentiment of DCASE 2018 ASC task 1 challenge has shown favor in the use of log mel-scaled spectrogram as the pre-processing technique and Convolution Neural Network (CNN) as the classifier and has shown promising results with top models achieving about 80% accuracy on DCASE 2018 Task 1. Hence, by following the leads from DCASE 2018, we created our "novel" model which address the important of both temporal information and frequency by using HPSS to best isolate harmonic and percussive component of an audio signal and further perform curated deep feature extraction using 'Time' specified CNN on the harmonic component and 'Frequency' specified CNN on the percussive component and finally combining this 2 features together using bilinear pooling. We present a 'twostream' architecture for sound classification. The remainder of this paper will be organized as such, Section II describe related works that inspire us to build our model, Section III elaborate on the feature extraction process, Section IV provide details on the model architecture, Section V explained about the experimental setup and dataset used. Section VI give a summary of the test results and Section VII conclude this paper.

II. MODEL DESIGN
In this paper, we present the use of harmonic percussive source separation (HPSS) [14] as the feature extraction technique and employing a CNN model on harmonic feature and a different CNN model on the percussive feature. Lastly, combining the two CNN models through bilinear pooling [7]. The choice of HPSS was inspired by the work of [5] and [6] who explored the used of HPSS for ASC tasks in DCASE 2018 and demonstrated the possibility of substituting HPSS over direct application of log-mel scale spectrogram.
Though the application of HPSS by [5] and [6] is the combining of harmonic feature and percussive feature as a 2channel feature input into a CNN model, this paper intuition is that harmonic feature and percussive feature should be processed by different CNN models with more affinity to the nature of the harmonic and percussive feature.
Hence, this paper explored the concept of a specified harmonic CNN (H-CNN) and percussive specified CNN (P-CNN) based on the idea driven by the work of [11], [14], [15] and [16]. These authors have a strong intuition that the current popular approached seen in DCASE 2018 of using a VGG [17] liked CNN as a classifier is not the best CNN model for sound classification tasks. Their reason behind this belief is that VGG liked CNN only employs a (3  3) convolution filter throughout its architecture, hence, with a small (3  3) local receptive field, it will impede the capability of the model in capturing the temporal relationship. (e.g., VGG which uses (3  3) filters will only capture the temporal relationship of 3 frames).
Of course, the simplest solution in solving this problem is by introducing a large convolution filter; unfortunately by doing so will also increase the computational complexity of the model. (e.g., 33 filters will have 9 weights and 77 filters will have 49 weights) Hence, these authors explored various filters size and concluded a "Time CNN" to have a long horizontal filter and "Frequency CNN" to have a long vertical filter. In conjunction with their findings, we proposed our H-CNN model to adopt a "Time CNN" filter strategy and P-CNN model to adopt a "Frequency CNN" filter strategy. Lastly, the effectiveness of combing these two models have been proven by [15], giving us the direction in combining these two models using bilinear pooling approach [7].

III. FEATURE EXTRACTION USING HPSS
This paper focus on exploiting the source separated from HPSS as input representation for the CNN model. HPSS was initially used on music classification where the idea is to separate harmonic sound instrument like, e.g. piano, wind instrument, and violin from a percussive sound instrument like e.g. drum and cymbal.
HPSS via median filtering is being selected as it is the most efficient approached [4]. The process of HPSS [4] in this setup can be broken down into four steps. Firstly, the Short-time Fourier transform (STFT) is applied to the raw audio signal to get a time-frequency representation. Secondly, median filtering with a kernel size of (1  31) is applied horizontally across the output from the first step to getting harmonic enhanced filtered spectrogram, and likewise, median filter with the similar configuration of (31  1) is also applied vertically to get percussive enhanced filtered spectrogram. Thirdly, with the output from the second step,, a masked binary is constructed, and the masked binary will be applied to the original STFT to get the harmonic component and percussive component. Lastly, both components are being processed into a log-mel spectrogram as shown in Fig. 1.
Unlike the usual HPSS performed on single audio to extract both harmonic and percussive features, this paper proposed a frequency biased HPSS to extract percussive feature and a time biased HPSS to extract harmonic feature. The difference between an HPSS for a percussive feature and the harmonic feature can be seen at the very first step, which is the STFT phase. By exploiting the weakness of STFT in regard to the tradeoff of bigger window size will result in a higher resolution in the frequency domain at the expense of losing some temporal details from the time domain. We proposed a huge window size of 1 sec to get the high resolution of frequency before performing median filtering and applying log-mel spectrogram with 256 filter banks. Thus, we determine this feature extraction processed as a frequency biased HPSS while time biased HPSS is treated in a general approach for the STFT with a window size of 42 ms and later applying log-mel spectrogram with 80 filter banks.
In this experiment, LibROSA [18] is being used as the package for HPSS and converting it to a log-mel spectrogram.

IV. CNN BILINEAR POOLING ARCHITECTURE
The proposed CNN architecture will consist of 2 CNN models, harmonic CNN, which will be described in A. moreover, percussive CNN which will be described in B.
These 2 CNN models will then be fused using bilinear pooling [7] and are named as Bilinear CNN (BCNN) by the author. Bilinear pooling is a technique that combines two CNN models by an outer matrix product on the local descriptors [19]: Equation (1) is a ( × ) vector, where = {ℎ 1 , ℎ 2 , … , ℎ ∈ ℝ } and likewise to = { 1 , 2 , … , ∈ ℝ } are each a set of local descriptors consisting of rows and columns, and S is the number of feature maps extracted through convolution. Local descriptors ℎ , are deep features extracted from H-CNN and P-CNN, respectively. Both CNN features are created from a standard convolution layers with a typical non-linear activation function Rectified Linear Unit (ReLU) [2], [1], [20] and usually follow by a widely used max pooling layer [20] to progressively reduce the spatial dimensionality in an effort to reduce computational power and also to reduce overfitting.
Hence, the output of H-CNN and P-CNN after the convolution layer is formulated as: Where I is either the input spectrogram or the feature maps with height of i and length of j and K is the filter which both will perform a matrix dot product while K convolute over I, denoted by ⨂, resulting in the output . K is a 2dimensional matrix ( , ) where n is the height of the filter and m is the length of the filter and = { 1 , 2 , … , } is a set of ( , )

Follow by batch normalization (BN) introduced in [21]
and applied on ASC task in [3] to tackle internal covariate shift cause by the increase variation between the weight parameters which will result in slower convergence time and possibility of vanishing or exploding gradient. Batch Normalization even out the output from (1) and later an activation function ReLU, ( ) = max (0, ) [1], [20] is applied, where x is all the features from the output of giving us the below equation: Lastly, max pooling which also have a filter KM with a dimension of (n, m) will convolute over the output while getting the maximum value of the features: = ⨂ = ∑ ∑ max ( ( , ), ( − , − )) (4) By combining the building blocks together, we get the output of h S and as depicted in the pseudo code in Table I  and Table II, respectively and the pseudo code of the proposed '2-stream' architecture in Table III. However, unlike a typical CNN architecture, instead of ending the convolution layers with a fully connected layer or max pooling, this paper proposed the use of bilinear pooling (1) with the intuition to capture pairwise correlations and features interaction, hence increasing the representative power of the deep features [7]. The proposed '2-stream' architecture approach in tackling ASC is illustrated in Fig. 2.

A. H-CNN (Harmonic CNN) Architecture
Similar to [11] and [15] our model will adopt the concept of Inception module [22] by employing a total of 10 horizontal filters of (1 × ), hence Q has a size of 10 and each set have a dimension of (1 × ) in (2), which will convolute around the input layer and later combining the features by depth-wise filter concatenation which is illustrated in Fig. 3 and depicted in Table I.   TABLE I. H-CNN PSEUDO CODE

Input:
Harmonic The intuition is to capture temporal relationship using various large horizontal filters from the harmonic source of the acoustic scene. The H-CNN model configuration is detailed in Table IV.

B. P-CNN (Percussive CNN) Architecture
Likewise, for P-CNN, the concept of Inception module is adopted, unlike H-CNN which specialized in capturing temporal relationship, P-CNN employs 10 vertical filters of ( × 1) , hence Q has a size of 10 and each set have a dimension of ( × 1) in (2), where n denotes the height of each filter. The idea here is to capture frequency which is the pitch aspect of the acoustic scene and the configuration of the CNN model is detailed in Table V. V. EXPERIMENTAL SETUP This paper uses dataset from DCASE 2019 Task 1a [23] which topic is on acoustic scene classification and researchers are given a dataset containing audio files to perform training and classification on 10 acoustic scenes. The 10 acoustic scenes being the airport, shopping mall, metro station, street pedestrian, public square, street traffic, tram, bus, metro and park.

A. Dataset
The development dataset for this task is the TAU Urban Acoustic Scenes 2019 Development dataset, is a collection of 10s audio clips which is recorded from 4 devices with a 48khz sampling rate and 24-bit resolution. The recording is conducted in 10 European cities and at 10 different location. A total of 40 hours of audio, with 14400 segments are collected in the development dataset and is evenly balanced between scene classes (144 segments per city per acoustic scene).
The trained model will then be evaluated on another 2 datasets, which is the TAU Urban Acoustic Scenes 2019 Leaderboard dataset and Evaluation dataset with Leaderboard dataset having a 3GB worth of 10s audio clips and Evaluation dataset containing 20 hours of audio recordings.

B. Training Setup and Evalution Setup
The training set consists of 9185 segments and 4185 segments for the testing set, which gives us the setup of splitting the development dataset by ~70% for training the model and ~30% for testing. For this challenge, the full 10s audio clip is being used, and no down sampling is being performed. Each audio clip is first, being converted from stereo into the mono channel and HPSS is being applied to split the audio into the harmonic wave and percussive wave.
Next, we preprocessed both waves into log-mel spectrogram with 50% overlapping, which will be the input representation of the network. The setup for the model training is as such; the model is trained with a mini-batch size of 16 samples over 200 epochs with patience of 15 epochs, meaning the model will stop training if validation loss did not improve after 15 epochs. Adam optimizer is being applied to the model with 0.001 learning rate, and learning rate will start to reduce if validation loss plateau for 5 epochs. The trained model is then being tested on the evaluation dataset.  Results show that our model performs slightly better than the DCASE 2019 baseline model and demonstrated that the proposed model is a feasible model for ASC tasks.
In a prior experiment, we have trained H-CNN model and P-CNN model separately with both models ending with a global max pooling follow by a softmax to determine the classification output. The accuracy performance achieved by both models is only between 57% to 59%. It means that both frequency and time features are equally important for our proposed model and experiment. Furthermore, this experiment result has proven that bilinear pooling is effective in combining frequency features and time features.

VII. CONCLUSION
To conclude, we have demonstrated that our proposed model has much potential and is now in its fetus stage. Hence, more work is required in fine-tuning the CNN architecture for both H-CNN and P-CNN model as currently, the number of filters is empirically decided and model enhancement technique such as dropout, data augmentation and hyperparameter training can be employed to further improve the model. Another improvement will be on exploring other bilinear pooling techniques. As the downside of employing bilinear pooling is its substantial computational cost which is also stated by [7] and [19]. In the feature extraction spectrum, we can investigate on other possible signal processing technique that establishes a more frequency domain bias and time domain bias.
Lastly, further acoustic classification test can be conducted to evaluate the robustness and capability of this model approach as this model should not be limited to only ASC tasks.