{'fraction_non_alphanumeric': 0.04424202992843201, 'fraction_numerical': 0.013012361743656473, 'mean_word_length': 4.228017132779038, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 34, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bridge Copula Model for Option Pricing

In this paper we present a new multi-asset pricing model, which is built upon newly developed families of solvable multi-parameter single-asset diffusions with a nonlinear smile-shaped volatility and an affine drift. Our multi-asset pricing model arises by employing copula methods. In particular, all discounted single-asset price processes are modeled as martingale diffusions under a risk-neutral measure. The price processes are so-called UOU diffusions and they are each generated by combining a variable (Ito) transformation with a measure change performed on an underlying Ornstein-Uhlenbeck (Gaussian) process. Consequently, we exploit the use of a normal bridge copula for coupling the single-asset dynamics while reducing the distribution of the multi-asset price process to a multivariate normal distribution. Such an approach allows us to simulate multidimensional price paths in a precise and fast manner and hence to price path-dependent financial derivatives such as Asian-style and Bermudan options using the Monte Carlo method. We also demonstrate how to successfully calibrate our multi-asset pricing model by fitting respective equity option and asset market prices to the single-asset models and their return correlations (i.e. the copula function) using the least-square and maximum-likelihood estimation methods.


Introduction
Many quantitative finance applications require a multi-asset pricing model with dependencies between the single-asset price components. Compared to the variety of univariate asset price models, the pool of multi-asset pricing models is not so extensive. Most of multivariate models are based on multidimensional geometric Brownian motion with the possible inclusion of a jump process. In this paper we develop and explore a new multi-asset arbitrage-free pricing model based on a special family of nonlinear diffusions. The development of efficient computational methods for pricing multi-asset equity derivatives under such a model and the calibration of the multi-asset model to both standard equity option data as well as historical equity prices are the objectives of the current paper.
Here, we specialize on option pricing applications under so-called UOU diffusion models, which are obtained by transforming an underlying Ornstein-Uhlenbeck diffusion process via the use of a diffusion canonical transformation method (see [2,3,5,8] and references therein). For all choices of model parameters, all discounted (single-asset) price processes UOU are conservative martingales under a riskneutral measure. Since the univariate diffusions are solvable, the single-asset risk-neutral transition probability density function is given in analytically closed form. Moreover, implied volatility surfaces for this highly nonlinear asset price model exhibit a wide range of pronounced smiles and skews of the type observed in the option markets. The main relevant features of the univariate UOU model are summarized in Section 1.
To construct a multivariate probability distribution, one can use a copula function that allows us to couple univariate distribution functions. Sampling from the obtained joint multivariate distribution function thereby reduces to sampling from the copula function and from the univariate distributions. Therefore, the copula method allows us to construct the joint distribution and density functions as well as to obtain an exact path sampling algorithm.
The main computational disadvantage of such an approach is the calculation of inverses of the distribution functions. This operation can be a rather time-consuming computational problem for a complicated multi-parameter distribution. Nevertheless, such a drawback can be significantly improved if the copula function and univariate distributions have a similar structure. As is shown in Subsection 1.3, the bridge probability density function (conditional on the values of the process at the endpoints of a time interval) of a UOU diffusion is reduced to a normal density. Hence it is natural to couple univariate UOU bridges using a Gaussian copula. Based on this idea, in Section 2, we construct a two-step algorithm for the exact path simulation of the multidimensional (nonlinear) UOU process in the risk-neutral measure. Firstly, we apply a usual copula method for sampling the multi-asset process at the terminal (maturity) time. Secondly, we use a bridge sampling along with a multivariate normal distribution to model the process at any intermediate time.
In Section 3, we demonstrate the calibration of the univariate and multivariate models to historical asset and equity option prices. The calibration process has two stages. First, we calibrate all univariate (marginal) asset price models independently of each other. Using the least-square method the models can be fitted to standard European option prices. Alternatively, the maximum likelihood estimation (MLE) allows the models to be fitted to historical asset prices. Second, we fit the copula function to historical observations. Since, our model assumes a normal copula, we need to find a best-fitted normal correlation matrix.
In Section 4, we give computational applications of the model to pricing multi-asset path-dependent Asian-style and Bermudan options. In pricing Bermudan options we use a regression-based Monte Carlo method.
In summary, the main results of our paper include: the construction a new family of multivariate models for which marginal processes are local volatility smile diffusions; the development of calibration schemes for the single-asset and multi-asset pricing UOU diffusion models based on the least-square and MLE methods; the construction of an exact multivariate path simulation method that can be used for Monte Carlo pricing of generally path-dependent European-style and American-style options.

Ornstein-Uhlenbeck Family of Univariate State-Dependent
Volatility Diffusion Models

Diffusion Canonical Transformation
The diffusion canonical transformation method, first presented in [2] and then further developed in [15,5,8], leads to various families of solvable one-dimensional diffusions with a nonlinear diffusion coefficient function and an affine drift. In this paper, we consider the UOU family, which is based on a regular Ornstein-Uhlenbeck process (X t ) t≥0 ∈ I ≡ R. The regular Ornstein-Uhlenbeck process is defined by the infinitesimal generator where λ > 0 and ν > 0 are constants. Both left and right boundaries l = −∞ and r = ∞ of the state space I are non-attracting natural for all choices of parameters. The transition probability density function (PDF) is where we define κ ≡ 2λ ν 2 > 0. Let ρ > 0 be a strictly positive constant. Then, two fundamental solutions of Lϕ(x) = sϕ(x), x ∈ I, are given by where υ ≡ ρ/λ > 0 and D −υ (·) is Whittaker's parabolic cylinder function (see [1]). The solutions ϕ + ρ (x) and ϕ − ρ (x) are linearly independent and respectively increasing and decreasing positive functions of x ∈ R.
We now construct another diffusion process (S t ) t≥0 ∈ D = R + = (0, ∞) by applying a diffusion canonical transformation to the underlying Ornstein-Uhlenbeck diffusion. This process obeys the stochastic differential equation (SDE) where r is a constant and σ(S) is a nonlinear diffusion coefficient function. Here (W t ) t≥0 is a standard Brownian motion in some appropriate probability measure P. The initial step of the transformation is to apply a Doob's h-transform or ρ-excessive transform (see [4]) to (X t ). The resulting diffusion process X (ρ) t t≥0 is defined by the following infinitesimal generator: where λ ρ (x) ≡ −λx + ν 2û ρ (x)/û ρ (x) for any ρ > 0. The strictly positive generating function is given byû The final step of the transformation is a change of variable (see [8] for a more general discussion). We define a new process S t = F(X (ρ) t ), t ≥ 0, that solves SDE (1.4) by finding a strictly monotonic map that solves L (ρ) F(x) = rF(x), for constant r. Then σ(F(x)) = ν|F (x)| or equivalently σ(s) = ν/|X (s)|, where X ≡ F −1 defines the unique inverse map. The transition PDF p S for the process (S t ) follows from that for the underlying process (X t ): x = X(s), x 0 = X(s 0 ). We now apply the above construction to a subfamily of diffusions with choice q 1 = 0, q 2 = 1, i.e. with generating functionû ρ (x) = ϕ − ρ (x). Letting r + ρ > 0 and ρ > 0, we specifically consider This function maps x ∈ R onto s ∈ (0, ∞) and is monotonically increasing, where s = F(x) has the unique inverse relation x = X(s). This transformation leads to a family of diffusions (S t ) that is referred to as the unbounded Ornstein-Uhlenbeck (UOU) family with the diffusion coefficient function given by where x = X(s). The volatility function σ in (1.9) depends on several adjustable positive parameters such as c, ν, λ, ρ, and more generally r = 0 also enters as an additional parameter within the volatility specification. Notice that for the driftless case with r = 0 formula (1.9) simplifies as follows: is the scale density for the X-diffusion.
The following is an important statement for the purposes of risk-neutral pricing.
. Consider a process (S t ) t≥0 ∈ R + of the UOU family solving the SDE (1.4) with the diffusion coefficient σ specified by (1.9) and having transition PDF p S specified by equations (1.7), (1.8), (1.9) and given by ] > 0 is the Wronskian 1 , and x = X(s), x 0 = X(s 0 ). Then (S t ) is a conservative process, i.e. P(S t ∈ R + ) = 1 for all t > 0. Moreover, the discounted asset price process (e −rt S t ) t≥0 is a true martingale.

Pricing Vanilla Options
According to Lemma 1.1 there exists an equivalent martingale measure under the UOU family with transition PDF (1.10) for any choice of the model parameters. Consider a standard European-style option defined by its payoff function Λ(S) at terminal price S = S T and maturity (expiration) time T > t 0 = 0. For example, a vanilla European call has the payoff function C E (S) = (S − K) + ≡ max{S − K, 0}, where K > 0 is a strike price. The valuation of a standard European option is given by the conditional expectation under a risk-neutral probability measure P ≡ Q: This is reduced to the valuation of a one-dimensional integral expressible as follows:  where x 0 = X(S 0 ).

The UOU Bridge
Let us consider a bridge process on [t 1 , t 2 ], 0 ≤ t 1 < t 2 , generated by a single asset price process S t with S t1 and S t2 fixed at s 1 and s 2 , respectively. The bridge density b S (t; s) ≡ b S (t 1 , t 2 , t; s 1 , s 2 , s) for S t = s, t 1 < t < t 2 , conditional on S t1 = s 1 and S t2 = s 2 , is given by is the bridge PDF of the underlying X-diffusion conditional on the endpoint values X t1 = x 1 and X t2 = x 2 . Since the underlying process (X t ) is a Gaussian process, the bridge process, generated by (S t ) t≥0 , with transition PDF (1.13) is just a nonlinear transformation of a Gaussian bridge, with respective path points mapped as X t = X(S t ).
Thus, the bridge PDF b S for (S t ) from the UOU family can be reduced to a Gaussian PDF for the underlying X-process. Indeed, from Eq. (1.2) we see that the bridge PDF b X of the Ornstein-Uhlenbeck diffusion is a normal density ϕ(x) = 1 √ 2πb e −(x−µ) 2 /2b 2 with mean µ and variance b 2 given by 2 Multivariate UOU-Diffusion Pricing Model

Coupling UOU Processes
Our goal is to construct a multi-asset price process (S t ) t≥0 with S t ≡ S 1 t , . . . , S n t , where each individual asset price process (S k t ) t≥0 , k = 1, 2, . . . , n, is a univariate UOU diffusion obeying (1.4) with common drift parameter r and diffusion function σ = σ k . Suppose that each of the n univariate processes is described by its own set of positive parameters ξ k = {λ k , ν k , c k , ρ k } , k = 1, 2, . . . , n. We denote by p k S and σ k the univariate risk-neutral transition PDF of the form (1.10) and the diffusion coefficient given by (1.9), respectively, which both correspond to the k-th asset price process. In the following, for each k = 1, 2, . . . , n,û k ρ k , F k , and X k will respectively denote the generating function, the mapping function, the inverse mapping function of the kth diffusion model. The transition PDFs p k X and p (ρ k ,k) X correspond to the underlying diffusion (X k t ) t≥0 and the transformed diffusion X . From the map, (F k ) (x) = σ k (F k (x))/ν k , we readily see that the respective correlations between the ratios of the infinitesimal increment and the diffusion coefficient remain the same after the change of variables: This relation shows how to couple UOU diffusions. First of all, we can directly couple the processes X (ρ k ,k) t , which, in turn, introduces correlations among the asset price processes (S k t ). In fact, the same method is used to couple geometric Brownian motions that are just exponentially transformed Brownian motions with drift. Hence, by using dependent Brownian motions, one introduces the correlations between the log-returns. In the case of more general diffusions, one may introduce correlations between the volatility-scaled returns. Alternatively, one may couple the underlying X-processes. For example, it is possible to use a multivariate Ornstein-Uhlenbeck process as an underlying vector process. Another approach (which is not discussed in this paper) is to extend the diffusion canonical transformation method directly to the case of multivariate diffusion processes.
Another idea, and the one that we implement in this paper, is the coupling of X (ρ) -bridges. As is shown in Subsection 1.3, the bridge PDF of a UOU X (ρ) -process is the same as the density function of the Ornstein-Uhlenbeck bridge, which is a multivariate normal PDF. The respective multivariate distribution function obtained with the use of a Gaussian copula is nothing more than a multivariate normal CDF. Therefore, the coupling of X (ρ) -bridges considerably simplifies the form of the joint path distribution function as well as the path simulation algorithm.

Copula Function
A copula C(u 1 , u 2 , . . . , u n ) is a multivariate cumulative distribution function (CDF) defined on the ndimensional unit cube [0, 1] n such that every marginal distribution is uniform on the unit interval [0, 1] (for a more detailed definition we refer to [16] . . , n. The well-known Sklar's theorem states that any n-dimensional joint distribution function Φ with continuous marginals Φ 1 , . . . , Φ n has a unique copula representation: In other words, for continuous multivariate distributions, the marginal distributions and the multivariate dependence structure can be separated. The multivariate density function f is then obtained by differentiating Eq. (2. 2): As a corollary of Sklar's theorem we have a class of copula functions constructed from continuous multivariate probability distributions as follows: The Gaussian (or normal) copula is one of the most important in financial applications. This copula is constructed from the multivariate normal distribution: where N R is the standard n-variate normal CDF with mean vector zero, unit variances, and correlation matrix R. Here, N −1 stands for the inverse of a standard univariate Gaussian CDF. Consider the problem of sampling from a multivariate distribution given by (2.2). The modelling algorithm consists of two steps. First, we simulate a uniformly distributed vector (U 1 , U 2 , . . . , U n ) ∈ [0, 1] n from the copula C. After that, we sample a vector (X 1 , X 2 , . . . , X n ) from the marginal distributions by evaluating the inverse CDFs: X k = (Φ k ) −1 (U k ), k = 1, 2, . . . , n. As is seen, the crucial part of this algorithm is the efficient computation of the inverse of a distribution function.

Multivariate Path Copula
The objective of this subsection is the derivation of a multivariate path distribution function in closed form. Such a function allows us to obtain an exact path simulation method and also to construct the joint path density function, which is used for calibrating the model with the maximum likelihood estimation method.
We employ the copula method to construct the joint distribution function of the multivariate process X is then obtained by applying the respective mapping function to each univariate X (ρ) -diffusion: Suppose that the process (X represent some arrangement of time points in T. The ordering of the time points is determined by the simulation method used. For the (forward) sequential method we assume that 0 In other words, we first obtain the value of the process at the terminal time T and then simulate a bridge path conditionally on the previously sampled value and the initial value. Notice that for the bridge sampling method, the ordering of intermediate time pointst 2 , . . . ,t N may vary, e.g. one can use a forward bridge sampling with 0 <t 2 <t 3 < · · · <t N <t 1 = T or a full bridge sampling method, where we successively halve the time interval or its largest segment to sample the process at the midpoint of the time segment conditionally on the previously sampled values at the endpoints of the current time (sub-)interval. The full bridge sampling may be applied along with the quasi-Monte Carlo method.
Let f k j (x) denote the PDF of X (ρ k ,k) tj conditional on the σ-algebra F k j−1 generated by the previously sampled j path points X x . Since all univariate processes are Markovian and are sampled backward in time, for each j = 2, 3, As is shown in Subsection 1.3, the PDF b k X is a normal density function with mean µ kj and variance b 2 kj , with values given by (1.14) where In X-space, the joint CDF Φ j and respective joint PDF f j of the n-dimensional point X From the Markov property of process (X In S-space, the joint path PDF f s and the joint PDFs f s j of St j , j = 1, 2 . . . , N , are respectively given by where s j = (s 1 j , . . . , s n j ) To sample from the distribution function in (2.6) we need a fast and accurate algorithm for inverting univariate CDFs. Usually, for non-standard distributions such as the PDF given by formula (1.10), such algorithms are quite computationally intensive. Therefore, the sequential approach becomes rather unfeasible for the path simulation of the multi-asset price UOU process (S t ).
As was mentioned above, we apply the Gaussian copula method to construct the underlying vector process X (ρ) t with transition PDFs of the form (1.6) instead of direct coupling of the S-space asset price processes. This is then followed by applying a nonlinear transformation to obtain the asset prices ). To minimize the number of numerical inversions of CDFs, we employ a bridge simulation approach, which is followed by the coupling of the bridge distributions. The application of the Gaussian copula to bridge PDFs from the UOU diffusion family leads to a multivariate normal distribution. Indeed, if each CDF Φ k (say, the CDF of the kth X (ρ) -bridge) is a normal CDF N x−µ k
Step 1. Apply the inverse mapping functions to obtain the initial values: Step 2. Sample the terminal-time value X N = X 1 N , . . . , X n N from the copula (2.2) by employing numerical inversion of the CDFs . . , n.
(i) Draw a normal vector (Z 1 , . . . , Z n ) from the n-variate normal distribution function N R with mean vector zero and correlation matrix R.
Step 4. Map the resulting values of the multivariate discrete-time process X (ρ) tj into the asset prices at each time point: As is seen from the algorithm, Step 1 involves the numerical inversion of n distribution functions. Since all parameters and initial spot prices are fixed, we can pre-compute and store the values of the inverse CDFs on a fine grid in (0, 1) and then apply a spline interpolation to sample X N (see [13]). Notice that the exact bridge simulation method presented above is faster than any approximation method such as the Euler scheme or the Milstein scheme (e.g., see [11]). First of all, our method has no limitations on time increments which can be very large. Therefore, if a path needs to be sampled only at a few time moments, no intermediate times are required to guarantee the convergence of sample paths. Second, the Euler approximation method (or any similar one) being applied to an SDE defined by (1.5) requires frequent computations of special functions that appear in the drift and diffusion coefficient functions. The bridge sampling method has only two computationally expensive steps, namely the sampling of terminal asset prices and the computation of mapping functions. As is mentioned above, the respective probability distributions can be tabulated to speed up the sampling. The mapping functions may be tabulated as well.
The approach presented here can be applied to other families of hypergeometric diffusions (see [7,8]). As is shown in [6], the probability distributions of the squared Bessel process and CIR process, which are the underlying diffusions for the so-called Bessel and confluent families of F -diffusions, respectively, are reduced to randomized gamma distributions. Such distributions are mixture probability distributions that can be obtained by allowing the shape parameter of the gamma distribution to be random. To couple F -diffusions, one can just couple the gamma distributions using a copula function. The use of a multivariate gamma distribution is more preferable, but all known examples of such a distribution admit only positive correlations.
3 Calibration of the Multivariate UOU Model

Univariate Case
It is very important from the practical point of view to develop a reliable and reasonably quick calibration scheme for the UOU diffusion family. Our objective is to obtain a calibration scheme that provides two levels of calibration: first, an initial full calibration of all parameters of the model and, second, a much faster recalibration that can be used as soon as new data have arrived. The second calibration scheme may be used throughout the day or even for longer periods, while the full calibration only need to be executed if markets move considerably.
To estimate a best-fitted parameter set ξ = {λ, ν, c, ρ} of the UOU model based on (observed) market option price data, the least squares method is employed. Suppose that an option with strike K i and maturity T i has an observed price O i , while the model produces a price of C i = C(K i , T i ; ξ) for the same option, where i = 1, 2, . . . , M . The goal of the calibration process is to minimize the least squares error for the M options considered: where w i is a weight that reflects the relative importance of reproducing the ith option price precisely. The suitable choice of the weight factors w i , i = 1, 2, . . . , M , is crucial for good calibration results. The confidence in individual data points is determined by the liquidity of the option. The weights can be evaluated from the bid-ask spreads: Alternatively, as it was suggested by [9], one may use the Black-Scholes (BS) "Vegas" evaluated at the implied volatilities of the market option prices to compute the weights: where ∂C BS /∂σ denotes the derivative of the BS option pricing formula with respect to the volatility σ, and σ BS In general, the calibration of a pricing model is an inverse problem, whose solution depends discontinuously on the data. To achieve uniqueness and stability of the solution, a penalty function is added to the least squares term: where the penalty function H is chosen such that the problem becomes well-posed.
As is examined in [9], the relative entropy method may be applied for solving ill-posed calibration problems. The relative entropy of a probability measure P on sample space Ω with respect to some primal measure P 0 is defined as follows: The regularization parameter α in (3.2) is used to adjust the trade-off between the accuracy of calibration and the numerical stability of results with respect to input option data. The right choice of α is based on the Morozov discrepancy principle [10], which is described by the following algorithm: 1. Compute parameters of ξ 0 of the primal measure P 0 by solving the nonlinear least squares problem (3.1) in low precision.

Numerical results for the univariate case
The data set used consists of 79 European call option prices with maturities ranging from less than one month up to 1.56 years. These market prices were obtained from Yahoo for IBM having the spot share value of 101.34 on July 7th, 2009. For the sake of simplicity, the risk-free interest is assumed to be constant and equal to r = 0.25%, and the dividend rate is set to zero. The calibration routine was developed using Matlab with the Optimization Toolbox, running on an Intel Core 2 CPU 2.14GHz with 2 GB of main memory.
To obtain the set of parameters for the primal probability measure, the UOU model is calibrated to the historical data from May 7th to July 7th, 2009. Using historical asset prices, S tj , j = 0, 1, . . . , N , 0 = t 0 < t 1 < · · · < t N , and the transition densities, we obtain the following (single-asset) log-likelihood function for this set of observations: Here, for simplicity, we assume the sequential simulation method. In case of a general sampling method, the log-likelihood function is given by where f j is defined analogously to f k j in Subsection 2.3 as if there was only one asset. In practice, the implementation of the calibration procedure is started with some initial values of parameters. The upper and lower bounds for the parameters should also be provided. Based on the empirical analysis, such bounds are obtained and are provided in Table 1.
The first step of the calibration procedure takes approximately 200 seconds to fit the model to 63 historical asset prices. The optimal values, that maximize the log-likelihood function (3.4), are ρ = 0.0357, υ = 0.0531, c = 118.2404, κ = 0.5951. This set of parameters defines the primal probability measure P 0 . The estimation of the regularization parameter α is based on the algorithm described above. The calculated value of α is 0.266.
The final step of the calibration process is the minimization of the nonlinear least squares function regularized by the relative entropy as is given in (3.2). The computation algorithm utilizes the Matlab function lsqnonlin with the exit tolerance set to 10 −6 . This function employs the Levenberg-Marquardt least-squares algorithm for estimating optimal parameters. The starting values and the limits for the parameters remain the same as given in Table 1. The computational algorithm takes approximately 400 seconds to fit the model to 79 option prices. The best-fitted parameters of the model are ρ = 0.0203, υ = 0.0013, c = 102.1384, κ = 0.6579. The objective function F α attains its minimum value of 1.58.
Notice that the discrepancy between the computed option prices and observed option prices may originate from different sources. First, the market data may contain errors or misleading information. For example, the values of illiquid options might be mispriced, or simple input errors may occur. Second, the calibration procedure estimates model parameters of an arbitrage-free model, while the market prices are not necessarily arbitrage-free. Hence, there is an inherent mismatch between the model prices and the market data. Notice that the use of time-dependent parameters may decrease the level and number of errors and make the calibration procedure maturity-wise. Another possible solution to improve the accuracy is to employ the calibration separately for out-of-the-money, at-the-money and in-the-money options.

Multivariate Case
Let us consider the multi-asset price processes (S t ) t≥0 modeled as described in Section 2, i.e. n UOU diffusions are coupled via the Gaussian copula function.
The calibration procedure can be split into two stages: (1) estimation of the parameters of the marginal (single-asset price) processes; (2) estimation of the correlation matrix R of the Gaussian copula. Such a calibration algorithm admits multiple variations. First, one may use the maximum likelihood estimation (MLE) to fit the marginal models to historical asset prices. Second, one may use the leastsquare method to fit the marginal models to historical derivative prices (say European options). For both approaches, the correlation matrix is then estimated by the MLE using historical asset prices. Alternatively, one may use only observed asset prices to estimate all parameter of the multivariate model simultaneously without splitting the calibration process. Notice also that the multivariate path distribution depends on the simulation method used. By using the sequential sampling or some version of the bridge sampling, one may obtain different models and, hence, obtain slightly different estimates of the model parameters.
be the n × (N + 1) matrix containing N + 1 independent historical prices for each of the n financial assets observed on a set of time points T = {t 0 , t 1 , . . . , t N }. Let (ξ, R) = (ξ 1 , . . . , ξ n , R) denote the set of parameters to be estimated. The historical observations in X (ρ) -space are obtained by applying the inverse map: X k tj = X k ( S k tj ; ξ k ). Suppose the joint path PDFs of the n-dimensional processes (X (ρ) t ) and (S t ) is constructed with the Gaussian copula as given by (2.6) and (2.7), respectively. The (n-asset) log-likelihood function is where φ R denotes the joint PDF of the n-variate normal distribution with mean vector zero, unit variances, and correlation matrix R; L 1 is the single-asset log-likelihood function given by (3.4) or (3.5), and L corr n denotes the log-likelihood function for the copula function. Recall that the expression in (3.6) is independent of the simulation method used. For the sequential and bridge methods, we provide below specific expressions of the log-likelihood function.
As is suggested by the structure of the log-likelihood function in (3.6), the calibration process can be split into two steps. First, the sets ξ k = {λ k , ν k , c k , ρ k } , k = 1, 2, . . . , n of parameters of the marginal distributions are estimated by employing the maximum likelihood estimation: As is seen, the parameters of the marginal distributions are estimated based on historical data. An alternative approach to computing the parameters is to fit asset price distributions to observed option prices as described in Subsection 3.1.
The last step is the estimation of the correlation matrix R for the given optimal univariate model parameters ξ ≡ { ξ 1 , . . . , ξ n }, ξ k ≡ λ k ,ν k ,ĉ k ,ρ k , estimated during the previous step.
First, we consider the sequential calibration method with the following log-likelihood function Sequential Calibration. For the sequential path generation method, the algorithm is as follows.
(ii) Compute vectors u j ≡ (u 1 j , . . . , u n j ) ∈ [0, 1] n , j = 1, . . . , N, by evaluating the integrals: (iii) Maximize the log-likelihood function with respect to R: The estimation of the log-likelihood function for the sequential calibration involves numerous estimations of the CDF for the UOU model. Since there is no simple-form solution for the CDF, the numerical integration of the probability density function should be performed regularly. By applying the bridge approach to the construction of the multivariate path distribution function, the number of integrals to be computed numerically on step (ii) reduces from n × N to n. This is due to the fact that for the bridge approach, the CDF Φ k j , j = 2, . . . , N , is Gaussian. Hence, for the bridge path generation method, the log-likelihood function can be simplified as follows: where mean µ kj and variance b 2 kj computed by formulae in (1.14); andt 1 = T .
Bridge Calibration. The following algorithm can be applied for the backward-in-time bridge path generation method.
(ii) Compute u N ≡ (u 1 N , . . . , u n N ), the values of normal CDFs corresponding to the terminal point of a path: (3.10) (iii) For each k = 1, . . . , n and j = 1, . . . , N − 1 calculate µ kj and b kj by using (1.14) with respective parameters λ = λ k , κ = κ k , ∆ 1 = t j − t 0 , and ∆ 2 = t j+1 − t j . Then, set (iv) Maximize the log-likelihood function with respect to R:   Table 2. Two approaches are then used for the evaluation of the optimal correlation matrix R. In the first approach, the correlation matrix is obtained by the pairwise calculation of the correlation coefficients. There are 4 2 correlation coefficients for 4-asset price processes to be estimated. In other words, instead of solving a 4-dimensional least-square problem with the log-likelihood function L corr n given by (3.8) or (3.9), we independently solve 4 2 2-dimensional problems with a 2-by-2 correlation matrix of the form 1 θ θ 1 and then compose another 4-by-4 correlation matrix using the estimations of θ.
However, the resulting matrix may violate the positive-definite property. To overcome this problem, a method suggested by [14] of finding the closest correlation matrix by the spectral decomposition is applied. The resulting matrices R are shown in Table 3. The computation time is 32.6 seconds for the bridge simulation and 47.9 seconds for the sequential simulation. In the second method, the correlation matrix as a whole is estimated. The computation of an optimal correlation matrix is performed in Matlab using the function fmincon, which allows us to find a minimum of a multivariate function with non-linear constraints. By adding nonlinear constraints, the algorithm works in the class of semi-positive matrices, which is absolutely necessary for the correct formulation of the correlation matrix. However, the candidate matrix, which minimizes the objective function in (3.6), may not have ones on the principal diagonal. To obtain a correct correlation matrix that is closest to the given one, the spectral decomposition method is applied again. The results are shown in Table 4 and Table 5.  Suppose the price process (S t ) is modelled at a discrete set of time points (written in an increasing order) T = {0 = t 0 , t 1 , t 2 , . . . , t N = T }. Hence, we construct a multivariate discrete-time n-dimensional price path (S i ) i=1,2,...,N , where S j = (S 1 j , S 2 j , . . . , S n j ). Let us consider two discrete-time monitored path-dependent securities, namely a Bermudan option and an Asian option. For an Asian-style option the payoff function Λ N (S 1 , . . . , S N ) is assumed to be a function of averages A k = A S k 1 , . . . , S k N of the asset prices, where, for example in the case of the arithmetic time-averaging, The Asian basket option with the payoff considered in Section 4.2, is an example of an arithmetic average option. The arbitrage-free value V = V (S 0 , T ) of a discrete-time monitored path-dependent option (without early exercise opportunities) can be represented as a mathematical expectation under a risk-neutral probability measure Q: Let V (t, S) denote the value of the option at time t given S t = S, assuming the option has not previously been exercised. We are interested in the present value V (t 0 , S 0 ). This value is determined through the backward-in-time recursion: where i = N − 1, . . . , 0. The risk-neutral expectation value function in the latter expression is called the continuation value function in state S at time t i−1 . In our computational tests, the univariate models correspond to the model parameter values used to plot the local volatility function with the thick line in Figure 1. So, for the UOU-family we simply set ρ k ≡ 0.02, κ k = 1, υ k = 0.5 and c k = 100, for all k. The interest rate is r = 0.05.

Pricing Asian Basket Options
In this example we use the Monte Carlo method for pricing the arithmetic Asian basket call option with payoff function (4.1) under the multivariate UOU model. The Monte Carlo simulations were run on the SHARCNET network (http://www.sharcnet.ca). Mostly we use Gulper-a 42-CPU cluster. The code was written in FORTRAN-90 using the MPI library.
First, we present numerical results for a bivariate (n = 2) UOU model. The correlation matrix used in the normal (Gaussian) copula has the usual form of R = 1 θ θ 1 . Table 6 provides the results of our numerical tests for the two-asset model. We used the following values of parameters: number of observations N = 100, terminal time T = 1.0, number of scenarios M = 10 000 000, and spot value S 0 = 100. The next example confirms that the computational complexity of the bridge copula method scales linearly with the number of assets. For all tests we use the crude Monte Carlo method on a single CPU. Clearly, the computational cost can be significantly reduced by using variance reduction methods and low-discrepancy point sets along with the parallel computing (see [7]). Another method of speeding up the computer code is to tabulate all distribution, generating and mapping functions that depend on the parabolic cylinder functions and other special functions. Table 7 provides the results of our numerical tests for increasing numbers of assets. We use the same values of parameters as those of the previous test but reduce the number of scenarios to M = 1 000 000. A 10-by-10 correlation matrix R given in Figure 3 is generated using the random Gram matrix method from [12]. For values of n < 10 we let R be the n × n-submatrix in the upper-left corner of the matrix given in Figure 3. Clearly, such a submatrix is again a correlation matrix of lower dimensionality. The idea of [12] is to generate independent pseudo-random vectors u 1 , . . . , u n distributed uniformly on the n-dimensional unit sphere and then to use the Gram matrix R = U U , where U ≡ (u 1 | · · · |u n ) has u k as kth column and U is the transpose of U . To create u k in R n we use a vector of independent standard normals z k ∼ N (0, I) normalized by the Euclidian norm: u k = z k / z k .

Pricing Bermudan Options
Regression-based methods are broadly applied to pricing multivariate American options (see [11] and references therein). The main idea consists in the use of regression to estimate continuation values from simulated paths. Each continuation value is approximated by a linear combination of some basis functions. The coefficients of such a representation are estimated by a regression method (typically least-squares). Table 8: Comparison of price estimates for Bermudan and European put options on the maximum of two assets modeled by the bridge copula multivariate UOU model with different values of the correlation coefficient θ. A sample standard error is given after ± sign.
To illustrate the regression algorithm above, we consider a Bermudan put option on the maximum of two underlying assets S 1 and S 2 modeled by the UOU bridge copula with interest rate r = 0.05, zero dividend yield, strike price K = 100, time to maturity T = 1.0, initial spot S 0 = 100, and N = 10 exercise times distributed evenly in the time interval [0, T ]. Here we use the same choice of the model parameters as in the previous numerical example. Clearly, this example can easily be extended to the case with three or more assets.
In the regression algorithm the basis functions are the power functions 1, S 1 , S 2 , S 2 1 , S 2 2 , S 1 S 2 , and payoff function Λ(S 1 , S 2 ) = (K − max(S 1 , S 2 )) + . We apply the regression method with M = 100 000 independent paths. To calculate standard error we replicate the random estimator 100 times. Table 8 shows numerical results for differing values of the correlation coefficient θ. For comparison's sake we also calculate the values of the European put option on the maximum of underlyings by the Monte Carlo method with 10 7 sample paths.

Conclusion
In this paper we have constructed a new multi-asset pricing model, whose marginal (single-asset price) processes are nonlinear (local volatility smile) regular diffusions from a recently developed UOU family of probability conserving martingale models [5,8]. The multivariate model is based on a bridge copula method, where a normal distribution function couples the underlying univariate Ornstein-Uhlenbeck bridges and consequently forms a multivariate asset price process with built-in correlations. Such an approach preserves the solvability of the model, hence the multivariate path density is available in closed form and can be used for the calibration of the model to market prices. Extra flexibility of the multivariate model is provided by different variations of the bridge sampling method. We are also able to sample multivariate paths from their exact distribution. Moreover, the proposed exact bridge simulation algorithm runs faster than any approximation scheme (e.g., the Euler method).
To illustrate the financial applications of our model, we succeeded in calibrating the model to singleasset equity option market prices as well as in calibrating the multi-asset price correlation matrix to historical asset prices. We also succeeded in pricing multi-asset Asian and Bermudan options by using a path integral Monte Carlo (MC) approach and an MC regression method, respectively. The preliminary calculations presented in this paper pave the way to further applications of derivatives pricing under such multi-asset state dependent diffusion models. The numerical efficiency of the path integral MC approach used in this paper can be further improved by employing quasi-MC methods. These methods have already been successfully used for another nonlinear model, the so-called Bessel-K model [5,8], and are also applicable to the UOU local volatility smile model.