Exact Simulation of IG-OU Processes

IG-OU processes are a subclass of the non-Gaussian processes of Ornstein–Uhlenbeck type, which are important models appearing in financial mathematics and elsewhere. The simulation of these processes is of interest for its applications in statistical inference. In this paper, a stochastic integral of Ornstein–Uhlenbeck type is represented to be the sum of two independent random variables—one has an inverse Gaussian distribution and the other has a compound Poisson distribution. And in distribution, the compound Poisson random variable is equal to a sum of Poisson-distributed number positive random variables, which are independent identically distributed and have a common specified density function. The exact simulation of the IG-OU processes, proceeding from time 0 and going in steps of time interval Δ, is achieved via the representation of the stochastic integral. Comparing to the approximate method, which is based on Rosinski’s infinite series representation of the same stochastic integral, by the quantile–quantile plots, the advantage of the exact simulation method is obvious. In addition, as an application, we provide an estimator of the intensity parameter of the IG-OU processes and validate its superiority to another estimator by our exact simulation method.      

Parallelization Method for a Continuous Property

An automated general purpose method is introduced for computing a rigorous estimate of a bounded region in ℝ ⁿ whose points satisfy a given property. The method is based on calculations conducted in interval arithmetic and the constructed approximation is built of rectangular boxes of variable sizes. An efficient strategy is proposed, which makes use of parallel computations on multiple machines and refines the estimate gradually. It is proved that under certain assumptions the result of computations converges to the exact result as the precision of calculations increases. The time complexity of the algorithm is analyzed, and the effectiveness of this approach is illustrated by constructing a lower bound on the set of parameters for which an overcompensatory nonlinear Leslie population model exhibits more than one attractor, which is of interest from the biological point of view. This paper is accompanied by efficient and flexible software written in C++ whose source code is freely available at .     

Balanced random interval arithmetic in market model estimation

The possibility of estimating bounds for the econometric likelihood function using balanced random interval arithmetic is experimentally investigated. The experiments on the likelihood function with data from housing starts have proved the assumption that distributions of centres and radii of evaluated balanced random intervals are normal. Balanced random interval arithmetic can therefore be used to estimate bounds for this function and global optimization algorithms based on this arithmetic are applicable to optimize it. The interval branch and bound algorithms with bounds calculated using standard and balanced random interval arithmetic were used to optimize the likelihood function. Results of the experiments show that when reliability is essential the algorithm with standard interval arithmetic should be used, but when speed of optimization is more important, the algorithm with balanced random interval arithmetic should be used which in this case finishes faster and provides good, although not always optimal, values.     

Error-free transformations of matrix multiplication by using fast routines of matrix multiplication and its applications

This paper is concerned with accurate matrix multiplication in floating-point arithmetic. Recently, an accurate summation algorithm was developed by Rump et al. (SIAM J Sci Comput 31(1):189–224, 2008). The key technique of their method is a fast error-free splitting of floating-point numbers. Using this technique, we first develop an error-free transformation of a product of two floating-point matrices into a sum of floating-point matrices. Next, we partially apply this error-free transformation and develop an algorithm which aims to output an accurate approximation of the matrix product. In addition, an a priori error estimate is given. It is a characteristic of the proposed method that in terms of computation as well as in terms of memory consumption, the dominant part of our algorithm is constituted by ordinary floating-point matrix multiplications. The routine for matrix multiplication is highly optimized using BLAS, so that our algorithms show a good computational performance. Although our algorithms require a significant amount of working memory, they are significantly faster than ‘gemmx’ in XBLAS when all sizes of matrices are large enough to realize nearly peak performance of ‘gemm’. Numerical examples illustrate the efficiency of the proposed method.     

Consistency tests in guaranteed simulation of nonlinear uncertain systems with application to an activated sludge process

In this paper, interval arithmetic simulation techniques are presented to determine guaranteed enclosures of the state variables of both continuous and discrete-time systems with uncertain but bounded parameters. In nonlinear uncertain systems axis-parallel interval boxes are mapped to complexly shaped regions in the state space that represent sets of possible combinations of state variables. The approximation of each region by a single interval box causes an accumulating overestimation from time-step to time-step, usually called the wrapping effect. The algorithm presented in this paper minimizes the wrapping effect by applying consistency techniques based on interval Newton methods. Subintervals that do not belong to the exact solution at a given time can be eliminated in order to give a tighter but still conservative approximation of the exact solution. Additionally, efficient splitting and merging strategies are employed to limit the number of subintervals. The proposed algorithm is applied to the simulation of an activated sludge process in biological wastewater treatment.     

Applications of interval arithmetic in solving polynomial equations by Wu''''s elimination method

Wu's elimination method is an important method for solving multivariate poly- nomial equations.In this paper,we apply interval arithmetic to Wu's method and convert the problem of solving polynomial equations into that of solving interval polynomial equa- tions.Parallel results such as zero-decomposition theorem are obtained for interval poly- nomial equations.The advantages of the new approach are two-folds:First,the problem of the numerical instability arisen from floating-point arithmetic is largely overcome.Second, the low efficiency of the algorithm caused by large intermediate coefficients introduced by exact compaction is dramatically improved.Some examples are provided to illustrate the effectiveness of the proposed algorithm.     

Provably recursive functions of constructive and relatively constructive theories

In this paper we prove conservation theorems for theories of classical first-order arithmetic over their intuitionistic version. We also prove generalized conservation results for intuitionistic theories when certain weak forms of the principle of excluded middle are added to them. Members of two families of subsystems of Heyting arithmetic and Buss-Harnik’s theories of intuitionistic bounded arithmetic are the intuitionistic theories we consider. For the first group, we use a method described by Leivant based on the negative translation combined with a variant of Friedman’s translation. For the second group, we use Avigad’s forcing method.     

An interval branch and bound method for global Robust optimization

In this paper, we design a Branch and Bound algorithm based on interval arithmetic to address nonconvex robust optimization problems. This algorithm provides the exact global solution of such difficult problems arising in many real life applications. A code was developed in MatLab and was used to solve some robust nonconvex problems with few variables. This first numerical study shows the interest of this approach providing the global solution of such difficult robust nonconvex optimization problems.

     

Fast evaluation of trigonometric polynomials from hyperbolic crosses

The discrete Fourier transform in d dimensions with equispaced knots in space and frequency domain can be computed by the fast Fourier transform (FFT) in arithmetic operations. In order to circumvent the ‘curse of dimensionality’ in multivariate approximation, interpolations on sparse grids were introduced. In particular, for frequencies chosen from an hyperbolic cross and spatial knots on a sparse grid fast Fourier transforms that need only arithmetic operations were developed. Recently, the FFT was generalised to nonequispaced spatial knots by the so-called NFFT. In this paper, we propose an algorithm for the fast Fourier transform on hyperbolic cross points for nonequispaced spatial knots in two and three dimensions. We call this algorithm sparse NFFT (SNFFT). Our new algorithm is based on the NFFT and an appropriate partitioning of the hyperbolic cross. Numerical examples confirm our theoretical results.     

An algorithm for automatically selecting a suitable verification method for linear systems

Several methods have been proposed to calculate a rigorous error bound of an approximate solution of a linear system by floating-point arithmetic. These methods are called ‘verification methods’. Applicable range of these methods are different. It depends mainly on the condition number and the dimension of the coefficient matrix whether such methods succeed to work or not. In general, however, the condition number is not known in advance. If the dimension or the condition number is large to some extent, then Oishi–Rump’s method, which is known as the fastest verification method for this purpose, may fail. There are more robust verification methods whose computational cost is larger than the Oishi–Rump’s one. It is not so efficient to apply such robust methods to well-conditioned problems. The aim of this paper is to choose a suitable verification method whose computational cost is minimum to succeed. First in this paper, four fast verification methods for linear systems are briefly reviewed. Next, a compromise method between Oishi–Rump’s and Ogita–Oishi’s one is developed. Then, an algorithm which automatically and efficiently chooses an appropriate verification method from five verification methods is proposed. The proposed algorithm does as much work as necessary to calculate error bounds of approximate solutions of linear systems. Finally, numerical results are presented.     

Roundoff error analysis of fast DCT algorithms in fixed point arithmetic

Discrete cosine transforms (DCT) are essential tools in numerical analysis and digital signal processing. Processors in digital signal processing often use fixed point arithmetic. In this paper, we consider the numerical stability of fast DCT algorithms in fixed point arithmetic. The fast DCT algorithms are based on known factorizations of the corresponding cosine matrices into products of sparse, orthogonal matrices of simple structure. These algorithms are completely recursive, are easy to implement and use only permutations, scaling, butterfly operations, and plane rotations/rotation-reflections. In comparison with other fast DCT algorithms, these algorithms have low arithmetic costs. Using von Neumann–Goldstine’s model of fixed point arithmetic, we present a detailed roundoff error analysis for fast DCT algorithms in fixed point arithmetic. Numerical tests demonstrate the performance of our results.      

Improved Gram–Schmidt Type Downdating Methods

The problem of deleting a row from a Q–R factorization (called downdating) using Gram–Schmidt orthogonalization is intimately connected to using classical iterative methods to solve a least squares problem with the orthogonal factor as the coefficient matrix. Past approaches to downdating have focused upon accurate computation of the residual of that least squares problem, then finding a unit vector in the direction of the residual that becomes a new column for the orthogonal factor. It is also important to compute the solution vector of the related least squares problem accurately, as that vector must be used in the downdating process to maintain good backward error in the new factorization. Using this observation, new algorithms are proposed. One of the new algorithms proposed is a modification of one due to Yoo and Park [BIT, 36:161–181, 1996]. That algorithm is shown to be a Gram–Schmidt procedure. Also presented are new results that bound the loss of orthogonality after downdating. An error analysis shows that the proposed algorithms’ behavior in floating point arithmetic is close to their behavior in exact arithmetic. Experiments show that the changes proposed in this paper can have a dramatic impact upon the accuracy of the downdated Q–R decomposition. AMS subject classification (2000) 65F20, 65F25     

Interval Newton Iteration in Multiple Precision for the Univariate Case

In this paper, interval arithmetic using an underlying multiple precision arithmetic is briefly presented. Then interval Newton iteration for solving nonlinear equations is introduced. A new Newton's algorithm based on multiple precision interval arithmetic is given, along with its properties: termination, arbitrary accuracy on the computed zeros, automatic and dynamic adaptation of the precision. Finally, some experiments illustrate the behaviour of this method.     

Computation of exact inertia and inclusions of eigenvalues (singular values) of tridiagonal (bidiagonal) matrices

This report may be considered as a non-trivial extension of an unpublished report by William Kahan (Accurate Eigenvalues of a symmetric tri-diagonal matrix, Technical Report CS 41, Computer Science Department, Stanford University, 1966). His interplay between matrix theory and computer arithmetic led to the development of algorithms for computing accurate eigenvalues and singular values. His report is generally considered as the precursor for the development of IEEE standard 754 for binary arithmetic. This standard has been universally adopted by virtually all PC, workstation and midrange hardware manufactures and tens of billions of such machines have been produced. Now we use the features in this standard to improve the original algorithm.In this paper, we describe an algorithm in floating-point arithmetic to compute the exact inertia of a real symmetric (shifted) tridiagonal matrix. The inertia, denoted by the integer triplet (π, ν, ζ), is defined as the number of positive, negative and zero eigenvalues of a real symmetric (or complex Hermitian) matrix and the adjective exact refers to the eigenvalues computed in exact arithmetic. This requires the floating-point computation of the diagonal matrix D of the LDL^t factorization of the shifted tridiagonal matrix T − τI with +∞ and −∞ rounding modes defined in IEEE 754 standard. We are not aware of any other algorithm which gives the exact answer to a numerical problem when implemented in floating-point arithmetic in standard working precisions. The guaranteed intervals for eigenvalues are obtained by bisection or multisection with this exact inertia information. Similarly, using the Golub-Kahan form, guaranteed intervals for singular values of bidiagonal matrices can be computed. The diameter of the eigenvalue (singular value) intervals depends on the number of shifts with inconsistent inertia in two rounding modes. Our algorithm not only guarantees the accuracy of the solutions but is also consistent across different IEEE 754 standard compliant architectures. The unprecedented accuracy provided by our algorithms could be also used to debug and validate standard floating-point algorithms for computation of eigenvalues (singular values). Accurate eigenvalues (singular values) are also required by certain algorithms to compute accurate eigenvectors (singular vectors).We demonstrate the accuracy of our algorithms by using standard matrix examples. For the Wilkinson matrix, the eigenvalues (in IEEE double precision) are very accurate with an (open) interval diameter of 6 ulps (units of the last place held of the mantissa) for one of the eigenvalues and lesser (down to 2 ulps) for others. These results are consistent across many architectures including Intel, AMD, SGI and DEC Alpha. However, by enabling IEEE double extended precision arithmetic in Intel/AMD 32-bit architectures at no extra computational cost, the (open) interval diameters were reduced to one ulp, which is the best possible solution for this problem. We have also computed the eigenvalues of a tridiagonal matrix which manifests in Gauss-Laguerre quadrature and the results are extremely good in double extended precision but less so in double precision. To demonstrate the accuracy of computed singular values, we have also computed the eigenvalues of the Kac₃₀ matrix, which is the Golub-Kahan form of a bidiagonal matrix. The tridiagonal matrix has known integer eigenvalues. The bidiagonal Cholesky factor of the Gauss-Laguerre tridiagonal is also included in the singular value study.     

On the proof of the reciprocity law for arithmetic Siegel modular functions

Earlier we obtained a new proof of Shimura’s reciprocity law for the special values of arithmetic Hilbert modular functions. In this note we show how from this result one may derive Shimura’s reciprocity law for special values of arithmetic Siegel modular functions. To achieve this we use Shimura’s classification of the special points of the Siegel space, Satake’s classification of the equivariant holomorphic imbeddings of Hilbert-Siegel modular spaces into a larger Siegel space, and, finally, a corrected version of some of Karel’s results giving an action of the Galois group Gal(Q_ab/Q) on arithmetic Siegel modular forms. Research supported in part by the NSF Grant No. DMS-8601130.     

Closed-form valuations of basket options using a multivariate normal inverse Gaussian model

This paper uses a multivariate normal inverse Gaussian model to develop closed-form pricing formulas for both geometric and arithmetic basket options. For geometric basket options, an exact analytical solution is possible; for arithmetic basket options, the formula is an approximation. The model is based on a jump-driven financial process, which is known empirically to be more realistic than a geometric Brownian motion. By comparing our results to Monte Carlo experiments, we confirm the internal consistency of our formulas. The “Greeks” can be derived from the closed-form formulas in a straightforward manner.     

An efficient algorithm for solving rank one perturbed linear Diophantine systems using Rosser��s approach

Recently, we described a generalization of Rosser’s algorithm for a single linear Diophantine equation to an algorithm for solving systems of linear Diophantine equations. Here, we make use of the new formulation to present a new algorithm for solving rank one perturbed linear Diophantine systems, based on using Rosser’s approach. Finally, we compare the efficiency and effectiveness of our proposed algorithm with the algorithm proposed by Amini and Mahdavi-Amiri (Optim Methods Softw 21:819–831, 2006).     

Adaptive and restarting techniques-based algorithms for circular packing problems

In this paper, we study the circular packing problem (CPP) which consists of packing a set of non-identical circles of known radii into the smallest circle with no overlap of any pair of circles. To solve CPP, we propose a three-phase approximate algorithm. During its first phase, the algorithm successively packs the ordered set of circles. It searches for each circle’s “best” position given the positions of the already packed circles where the best position minimizes the radius of the current containing circle. During its second phase, the algorithm tries to reduce the radius of the containing circle by applying (i) an intensified search, based on a reduction search interval, and (ii) a diversified search, based on the application of a number of layout techniques. Finally, during its third phase, the algorithm introduces a restarting procedure that explores the neighborhood of the current solution in search for a better ordering of the circles. The performance of the proposed algorithm is evaluated on several problem instances taken from the literature.     

Computational existence proofs for spherical <Emphasis Type="Italic">t</Emphasis>-designs

Spherical t-designs provide quadrature rules for the sphere which are exact for polynomials up to degree t. In this paper, we propose a computational algorithm based on interval arithmetic which, for given t, upon successful completion will have proved the existence of a t-design with (t + 1)² nodes on the unit sphere S² ì \mathbbR³{S^2 \subset \mathbb{R}^3} and will have computed narrow interval enclosures which are known to contain these nodes with mathematical certainty. Since there is no theoretical result which proves the existence of a t-design with (t + 1)² nodes for arbitrary t, our method contributes to the theory because it was tested successfully for t = 1, 2, . . . , 100. The t-design is usually not unique; our method aims at finding a well-conditioned one. The method relies on computing an interval enclosure for the zero of a highly nonlinear system of dimension (t + 1)². We therefore develop several special approaches which allow us to use interval arithmetic efficiently in this particular situation. The computations were all done using the MATLAB toolbox INTLAB.     

The theorems of Loewner and Pick

The study of the exact interpolation of quadratic norms in vector spaces depends in an essential way on the theory of monotone matrix functions developed by Loewner in 1934 [4]. This theory, in its turn, depends on Loewner’s solution of a problem of interpolation by rational functions of a certain class. The discussion of this latter problem is necessarily complicated, and Loewner’s text does not lend itself to ready reference. It has therefore seemed worthwhile to recast a portion of Loewner’s results in a form more suited to the applications we have in view. Our work, however, is not wholly derivative; none of our theorems are explicitly stated by Loewner and our arguments, which are of a more geometric character, are essentially different. The knowledgeable reader will note that our hypotheses are slightly stronger than Loewner’s and that our results are therefore also stronger. For the applications which we have in mind, Theorem III is the most important result; the proof of this theorem depends on all of the previously developed theory.