Dalal Sukkari

Hatem Ltaief

Extreme David Keyes

Mathieu Faverge
email: mathieu.faverge@inria.fr

Asynchronous Task-Based Polar Decomposition on Manycore Architectures

Keywords: Polar decomposition, Asynchronous execution, Dynamic runtime system, Fine-grained execution, Directed acyclic graph, High performance computing 0.1

This paper introduces the first asynchronous, taskbased implementation of the polar decomposition on manycore architectures. Based on a new formulation of the iterative QR dynamically-weighted Halley algorithm (QDWH) for the calculation of the polar decomposition, the proposed implementation replaces the original and hostile LU factorization for the condition number estimator by the more adequate QR factorization to enable software portability across various architectures. Relying on fine-grained computations, the novel task-based implementation is also capable of taking advantage of the identity structure of the matrix involved during the QDWH iterations, which decreases the overall algorithmic complexity. Furthermore, the artifactual synchronization points have been severely weakened compared to previous implementations, unveiling look-ahead opportunities for better hardware occupancy. The overall QDWH-based polar decomposition can then be represented as a directed acyclic graph (DAG), where nodes represent computational tasks and edges define the inter-task data dependencies. The StarPU dynamic runtime system is employed to traverse the DAG, to track the various data dependencies and to asynchronously schedule the computational tasks on the underlying hardware resources, resulting in an out-of-order task scheduling. Benchmarking experiments show significant improvements against existing stateof-the-art high performance implementations (i.e., Intel MKL and Elemental) for the polar decomposition on latest shared-memory vendors' systems (i.e., Intel Haswell/Broadwell/Knights Landing, NVIDIA K80/P100 GPUs and IBM Power8), while maintaining high numerical accuracy.

I. INTRODUCTION

Today's most powerful supercomputers are composed of fat computational nodes over-provisioned by floating-point units [START_REF]The Top500 List[END_REF], which may distort the balance of characteristics systems with respect to other hardware resources (the Kiviat diagram), such as memory per core, aggregated bandwidth, I/O nodes, interconnect, etc. Although real scientific applications are often memory-bound with low arithmetic intensity kernels, and therefore limited by the bus bandwidth, we revisit the polar decomposition, an important dense linear algebra (DLA) algorithm, which can make an effective use of the predominant floating-point units provided by the current stateof-the-art hardware vendor chips (for instance, Intel Knights Landing and NVIDIA Pascal P100). Based on the QR-based dynamically weighted Halley (QDWH) iteration introduced by Nakatsukasa et. al [START_REF] Nakatsukasa | Optimizing Halley's Iteration for Computing the Matrix Polar Decomposition[END_REF], the polar decomposition is a key algorithm for various scientific applications, e.g., in continuum mechanics to decompose stress tensors and to simulate the deformation of an object, in aerospace computations [START_REF] Meyer | Practical Comparison of Iterative Matrix Orthogonalization Algorithms[END_REF] during strapdown inertial navigation and other aerospace systems to describe the rotation of one coordinate system relative to a reference coordinate system, and in chemistry [START_REF] Goldstein | Linear algebra and quantum chemistry[END_REF] to help the understanding of properties in terms of electron pair (chemical bond) transferability, etc. Further applications are also reported by Higham in [START_REF] Higham | Computing the Polar Decompositionwith Applications[END_REF].

This paper describes the first asynchronous, task-based implementation of the QDWH-based polar decomposition on manycore architectures. The standard algorithm requires up to six iterations to converge and to calculate the polar factor, depending on the condition number of the input matrix, involving O(n 3 ) matrix operations at each operation. Its algorithmic complexity may, therefore, be prohibitive. Nevertheless, this challenge can be compensated for the high level of concurrency exposed at each iteration [START_REF] Sukkari | A High Performance QDWH-SVD Solver Using Hardware Accelerators[END_REF], [START_REF]High Performance Polar Decomposition on Distributed Memory Systems[END_REF].

This paper proposes to considerably improve previous works [START_REF] Nakatsukasa | Optimizing Halley's Iteration for Computing the Matrix Polar Decomposition[END_REF], [START_REF] Sukkari | A High Performance QDWH-SVD Solver Using Hardware Accelerators[END_REF], [START_REF]High Performance Polar Decomposition on Distributed Memory Systems[END_REF] from two distinct algorithmic and implementation perspectives. The former consists in replacing the hostile LU-based matrix condition number estimation by an adequate QR-based implementation for broader code portability across vendor architectures. The latter has twofold aspects: [START_REF]The Top500 List[END_REF] it permits to take advantage and exploit the structure of the identity matrix involved at each QR-based QDWH iteration, which significantly reduces the algorithmic complexity, thanks to fine-grained computations associated with a dynamic asynchronous execution, and (2) the artifactual synchronization points are severely weakened, unveiling lookahead opportunities for better hardware occupancy. The overall QDWH-based polar decomposition can then be represented as a directed acyclic graph (DAG), where nodes represent computational tasks and edges define the inter-task data dependencies. We employ the StarPU dynamic runtime system to unroll the DAG, to track the various data dependencies and to asynchronously schedule the computational tasks on the underlying hardware resources, resulting in an out-oforder task scheduling. StarPU increases user-productivity by establishing a separation of concerns consisting in hiding the hardware complexity from library developers. This enables end-users to target various hardware architectures with a single source code. Extensive benchmarking experiments show significant improvements against existing state-of-the-art high performance implementations (i.e., MKL and Elemental) for the polar decomposition on latest shared-memory systems (i.e., Intel Haswell/Broadwell/Knights Landing, NVIDIA K80/P100 GPUs and IBM Power8), while maintaining high numerical accuracy for well and ill conditioned matrices.

The remainder of the paper is organized as follows. Section II presents related work. Section III highlights our research contributions. Section IV briefly recalls the polar decomposition and its main computational phases. Section V describes the fundamental design of current state-of-the-art DLA software libraries, as implemented in LAPACK [START_REF] Anderson | LAPACK User's Guide[END_REF], MAGMA and PLASMA [START_REF] Agullo | Numerical Linear Algebra on Emerging Architectures: The PLASMA and MAGMA projects[END_REF]. The implementation details of the high performance task-based asynchronous QDWH are given in Section VI. Section VII provides new upper-bound for the QDWH algorithmic complexity. Numerical accuracy, implementations assessments and performance comparisons with existing state-of-the-art DLA softwares are given in Section VIII and we conclude in Section IX.

II. RELATED WORK

The polar decomposition algorithm has been well studied in the last three decades in terms of complexity and numerical robustness/accuracy [START_REF] Gander | On Halley's iteration method[END_REF]- [START_REF] Laszkiewicz | Approximation of matrices and a family of gander methods for polar decomposition[END_REF]. It consists in decomposing a dense matrix A = U p H, where U p is the orthogonal polar factor and H is the positive semi-definite Hermitian polar factor. Initially designed with Newton's method based on an explicit matrix inversion calculation, numerical instability was reported, especially in presence of ill-conditioned matrices. An algorithm based on Halley's iteration was then introduced with asymptotically cubic rate of convergence in obtaining the final polar factor. To solve the numerical accuracy issues due to the matrix inversion computation, an inverse-free QR-based dynamically weighted Halley (QDWH) has finally been proposed by Nakatsukasa et. al [START_REF] Nakatsukasa | Optimizing Halley's Iteration for Computing the Matrix Polar Decomposition[END_REF]. However, the polar decomposition has not been implemented in a high performance computing system's environment, most probably due to its excessive algorithmic complexity, which does not reflect a practical assessment of the method. More recently, Nakatsukasa and Higham [START_REF] Nakatsukasa | Stable and Efficient Spectral Divide and Conquer Algorithms for the Symmetric Eigenvalue Decomposition and the SVD[END_REF] have shown that QDWH can be used as a building block for the dense symmetric eigensolver and singular value decomposition [START_REF] Golub | ser. John Hopkins Studies in the Mathematical Sciences[END_REF], [START_REF] Trefethen | Numerical Linear Algebra[END_REF], which has brought to the fore further research directions. Indeed, previous works from the authors have implemented QDWH-based singular value decomposition on hardware accelerators [START_REF] Sukkari | A High Performance QDWH-SVD Solver Using Hardware Accelerators[END_REF] and distributed-memory systems [START_REF]High Performance Polar Decomposition on Distributed Memory Systems[END_REF], where the calculation of the polar factor is the mosttime consuming phase. The aforementioned implementations have somewhat demonstrated limited performance scalability on multiple GPUs and large clusters. This is mostly due to the low hardware resource occupancy achieved by the inherent bulk synchronous programming model (BSP), which both implementations rely on for parallel performance. By the same token, it is also noteworthy to mention that the high performance software library Elemental [START_REF] Poulson | Elemental: A new framework for distributed memory dense matrix computations[END_REF] provides a QDWH implementation for distributed-memory systems.

Last but not least, the polar decomposition can alternatively be computed through an SVD as follows:

A = U ΣV = U V V ΣV = U p V ΣV = U p H.
This strategy has shown some performance scalability issues, due to the slow convergence of the QR algorithm on the condensed bidiagonal form [START_REF]High Performance Polar Decomposition on Distributed Memory Systems[END_REF].

III. CONTRIBUTIONS

The following contributions represent the crux of the paper: (1) improve the standard QDWH algorithm by replacing the LU-based condition estimator with QR, without increasing the overall algorithmic complexity, while enabling software portability across hardware architectures, (2) develop the first task-based QDWH implementation based on fine-grained computations, which enables to exploit the identity data structure during the QDWH iterations, reducing up to 20% the overall complexity, (3) rely on a dynamic runtime system (i.e., StarPU) to asynchronously schedule the computational tasks among available processing units in order to improve hardware occupancy, (4) provide a comprehensive performance assessment and comparisons on a myriad of high-end architectures.

IV. THE POLAR DECOMPOSITION

The paper focuses on the inverse-free QDWH-based iterative procedure to calculate the polar decomposition [START_REF] Nakatsukasa | Optimizing Halley's Iteration for Computing the Matrix Polar Decomposition[END_REF], [START_REF] Nakatsukasa | Stable and Efficient Spectral Divide and Conquer Algorithms for the Symmetric Eigenvalue Decomposition and the SVD[END_REF] of a matrix A ∈ R m×n (m ≥ n), such that A = U p H. To ensure the paper is self-contained, we briefly recall the convergent sequence as follows, with A the initial matrix:

U 0 = A/α, α = A 2 U k+1 = b k c k U k + 1 √ c k a k -b k c k Q 1 Q 2 , k ≥ 0 (1)
where U p = lim k→∞ (U k ), and

Q 1 Q 2 R = √ c k U k I .
H can then be found with the two steps formula:

H = U p * A, H = 1 2 (H + H ) (2)
The main goal consists in calculating the optimal parameters (a k , b k , c k ) so that cubical convergence can be attained during the QDWH iteration. The expression of the parameters (a k , b k , c k ) can be written as follows:

a k = h(l k ), b k = (a k -1) 2 /4, c k = a k + b k -1, l 0 = β α , l k = l k-1 (a k-1 + b k-1 l 2 k-1 ) 1 + c k-1 l 2 k-1 , k >= 1, h(l) = √ 1 + d + 1 2 8 -4d + 8(2 -l 2 ) l 2 √ 1 + d , d = 3 4(1 -l 2 ) l 4 , (3)
with β = 1/ A -1 2 . For ill-conditioned matrices, the number of iterations k can vary up to six. We refer to [START_REF] Nakatsukasa | Optimizing Halley's Iteration for Computing the Matrix Polar Decomposition[END_REF] for further details on the theoretical proof. When U k becomes well-conditioned, it is possible to replace Equation 1 with a Cholesky-based implementation as follows:

U k+1 = b k c k U k + a k - b k c k (U k W -1 k )W - k , W k = chol(Z k ), Z k = I + c k U k U k . (4)
This algorithmic switch at runtime allows to further speed up the overall computation, thanks to a lower algorithmic complexity, while still maintaining numerical stability. In practice, this transition is monitored by setting a threshold for c k . Once convergence is reached, the polar factor is U p = U k and the positive semi-definite Hermitian polar factor corresponds to H = U p A. All in all, the number of floating-point operations depends on the number of iterations required to converge, which is dictated by the condition number of the original matrix problem. Typically, for illconditioned matrices, QDWH will perform three QR-based QDWH iterations (Equation 1), followed by three Choleskybased QDWH iterations (Equation 4), besides executing other compute-intensive Level 3 BLAS operations, i.e., triangular solves, applications of Householder reflectors, matrix-matrix multiplications etc.

V. HIGH PERFORMANCE DLA SOFTWARE LIBRARIES

As discussed in Section IV, although the QDWH-based polar decomposition is a challenging and complex algorithm, it relies on conventional dense linear algebra operations, e.g., QR/Cholesky-based linear solvers, which are well-supported by several open-source and vendor-optimized numerical libraries. These libraries can be differentiated into the two following algorithmic categories.

A. Block Algorithms

Block algorithms rely on successive panel and update sequences to perform matrix computations. The panel phase is memory-bound and does not benefit from thread parallelism, while the phase of the trailing submatrix update is highly parallel, in which computations are applied by means of multithreaded Level 3 BLAS kernel executions. These sequences are characteristic of the fork-join paradigm, alternating sequential and parallel computational phases, and therefore, suffer from performance losses due to low hardware occupancy engendered by unnecessary in-between synchronization points. In fact, this bulk synchronous programming model corresponds to the backbone of many open-source and vendors' state-of-the-art numerical libraries such as LAPACK [START_REF] Anderson | LAPACK User's Guide[END_REF], MAGMA [START_REF] Magma | Matrix Algebra on GPU and Multicore Architectures[END_REF] and ScaLAPACK [START_REF] Blackford | ScaLAPACK Users' Guide[END_REF] for shared-memory, accelerator-based and distributed-memory systems, respectively. As highlighted in the exascale software roadmap [START_REF] Dongarra | The International Exascale Software Project Roadmap[END_REF], BSP models may need to be reconsidered, especially in presence of millions of cores, which already constitute today's supercomputers [START_REF]The Top500 List[END_REF].

B. Tile Algorithms

To answer this call for action and provide a solution for the challenge brought by the manycore era, the DLA community has initiated a decade ago a profound redesign of matrix computation algorithms in order to benefit from the high level of concurrency. This translated into breaking down the dense matrix data structure into tiles following a tile data layout as opposed to the standard column-major format, which is the standard for block algorithms. The various matrix operations can then be represented as a directed acyclic graph (DAG), where nodes represent sequential computational tasks and edges define the inter-task data dependencies. The resulting fine-grained computations permit to weaken the artifactual synchronization points by bringing to the fore lookahead opportunities, which in return, can be exploited by dynamic runtime systems in keeping threads in a busy state throughout the entire execution. The performance gain of tile versus block algorithms have been thoroughly addressed in the literature [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF]- [START_REF] Agullo | Comparative study of one-sided factorizations with multiple software packages on multicore hardware[END_REF], in the context of PLASMA [START_REF]PLASMA Users' Guide, Parallel Linear Algebra Software for Multicore Architectures[END_REF] and FLAME [START_REF]The FLAME project[END_REF] numerical software libraries, using QUARK [START_REF] Yarkhan | QUARK Users' Guide: QUeueing And Runtime for Kernels[END_REF] and SuperMatrix [START_REF] Chan | Supermatrix out-of-order scheduling of matrix operations for smp and multi-core architectures[END_REF] runtime systems, respectively. More recently, in a community effort to enhance user productivity by abstracting the hardware complexity, the Chameleon library [START_REF]The Chameleon project[END_REF] has been developed to target multiple hardware architectures with a single source code. This is achieved by standardizing existing dynamic runtime system APIs (e.g., OpenMP [START_REF] Openmp | Openmp application program interface version 4[END_REF], OmpSs [START_REF] Duran | A Proposal to Extend the OpenMP Tasking Model with Dependent Tasks[END_REF], QUARK [START_REF] Yarkhan | QUARK Users' Guide: QUeueing And Runtime for Kernels[END_REF], StarPU [START_REF] Augonnet | StarPU: A unified platform for task scheduling on heterogeneous multicore architectures[END_REF], PaRSEC [START_REF] Bosilca | DAGuE: A generic distributed DAG engine for High Performance Computing[END_REF]) through a thin layer of abstraction, making the user developer experience oblivious to the underneath runtime system and its corresponding hardware deployment. For instance, this oblivious software infrastructure has been already used in the context of computational astronomy using StarPU [START_REF] Charara | Pipelining Computational Stages of the Tomographic Reconstructor for Multi-Object Adaptive Optics on a Multi-GPU System[END_REF], and more recently with OmpSs [START_REF] Ltaief | Adaptive Optics Simulation for the World's Largest Telescope on Multicore Architectures with Multiple GPUs[END_REF].

The QDWH-based polar decomposition resurfaces during a possible period of convergence, where hardware/software codesign plays now a major role in designing future systems and numerical libraries for exascale.

VI. HIGH PERFORMANCE IMPLEMENTATIONS

In this section, we describe the task-based implementation of the QDWH algorithm and the novel optimizations introduced to increase hardware occupancy and overall performance, in the context of the Chameleon library [START_REF]The Chameleon project[END_REF].

A. The StarPU Runtime System

StarPU is the de facto dynamic runtime system for Chameleon. Although Chameleon supports other runtimes (e.g., PaRSEC [START_REF] Bosilca | DAGuE: A generic distributed DAG engine for High Performance Computing[END_REF], QUARK [START_REF] Yarkhan | QUARK Users' Guide: QUeueing And Runtime for Kernels[END_REF] and recently OmpSS [START_REF] Duran | A Proposal to Extend the OpenMP Tasking Model with Dependent Tasks[END_REF]), we decided to solely rely on the StarPU [START_REF] Augonnet | StarPU: A unified platform for task scheduling on heterogeneous multicore architectures[END_REF] runtime system to implement this algorithm, since it is probably one of the most mature runtime systems when it comes to supporting various hardware architectures. StarPU is a runtime, which deals with the execution of generic task graphs. This task graph is given through the sequential task flow (STF) programming model where tasks are inserted to the runtime in a sequential manner with additional hints on the data usage: read, write, read-write. Then, the runtime automatically infers data dependencies from those hints while unrolling the sequential flow of task submissions, and is in charge of scheduling the tasks while enforcing those dependencies. The PaRSEC [START_REF] Bosilca | DAGuE: A generic distributed DAG engine for High Performance Computing[END_REF] runtime is also used in dense linear algebra libraries but relies on a parameterized task graph (PTG) representation of the application. This model allows for compact and problem size independent representations of the graph to execute but is usually less intuitive to the non familiarized user.

One of the main advantages of using the task-based implementation is to become oblivious of the targeted architectures. This improves the user productivity, and it is even more realistic for runtimes as StarPU, which are able to transparently handle single heterogeneous nodes, and eventually multiple heterogeneous nodes in case the StarPU-MPI [START_REF] Agullo | Achieving High Performance on Supercomputers with a Sequential Task-based Programming Model[END_REF] extension is used. To enable such portability, StarPU tasks are associated to codelets which groups under the same name multiple implementations of the same task: CPU, CUDA, OpenCL, OpenMP, . . . At runtime, StarPU will automatically decide which implementation of the task is better suited to achieve the highest performance based on cost models. These cost models are automatically generated by StarPU when executing the application and kept for subsequent executions. These models are especially important to the Heterogeneous First Time [START_REF] Topcuoglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF] (HeFT) scheduling strategy used by StarPU, when accelerators are involved in the computations. Another benefit from using such programming models is the capabilities offered to the programmer to submit simultaneously independent steps of an application to raise the resources occupancy, and add a single synchronization point when all steps are performed. The MORSE_xxxx_Tile_Async interface of the Chameleon library offers this capability to interleave multiple dense linear algebra operations when it is possible. Conversely, the non-asynchronous interface, MORSE_xxxx_Tile, enforces a synchronization call at the end of the function to wait for the end of all submitted tasks.

B. Task-Based QDWH Polar Decomposition Pseudo-Code

Algorithm 1 presents the pseudo-code of the task-based QDWH implementation on top of the Chameleon library. It is decomposed in three main code sections. The first one from row 1 to 6 evaluates the two-norm of the input matrix A that is required to start the iterative process. The specificity of this code section is the 2-norm estimator genm2 that we introduced in the Chameleon library through an iterative computation in which we tried to minimized the number of synchronizations. The second section of the algorithm computes the initial condition number l 0 from row 7 to 19. The classical way consists in computing an LU factorization of the matrix A, and its one-norm. Then, it is possible to compute an estimator of the condition number with those two information (dgecon). The main challenge here resides in the LU factorization with partial pivoting, which is difficult to implement using task-based programming model, due to the global synchronization points needed during the panel factorization while looking for pivot candidates and the resulting row swapping step. Some solutions have been proposed on sharedmemory systems [START_REF] Dongarra | Exploiting finegrain parallelism in recursive lu factorization[END_REF] but there are no existing solutions that are oblivious of heterogeneous architectures. We thus propose

L2 = Li 2 , dd = 3 (4(1 -L2)/L2 2 )
25:

sqd = √ 1 + dd 26: a1 = sqd + 8 -4 × dd + 8(2 -L2)/(L2 × sqd))/2
27:

a = real(a1), b = (a -1) 2 /4, c = a + b - 1
28:

Li = Li(a + b × L2)/(1 + cL2)

29:

dlacpy Async( U, U1 ) Backup U k-1 30: 31:

/* Compute U k from U k-1 */

32:

if c > 100 then 33:

C = C1 C2 = √ cU k-1 I 34: dgeqrf Async( C ) C = QR = Q1 Q2 R 35: dorgqr Async( C ) C = Q = Q1 Q2 36: dgemm Async( Q1, Q 2 , U ) U k = 1 √ c a -b c Q1Q 2 + b c U k-1

37:

else 38:

dlaset Async( Z, 0., 1. ) Z = I 39:

dgemm Async( U , U , Z) Z k = I -cU k-1 U-1k 40: dgeadd Async( U , B ) B = U k-1 41: dposv Async( Z, B ) Solve Z k x = U k-1 42: dgeadd Async( B, U ) U k = b c U k-1 + a -b c (U k-1 W -1 k-1 )W k-1

43:

end if

44:

dgeadd Async( U, U1 )

U k -U k-1 45: dlange Async( U1, conv ) conv = U k -U k-1 F

46:

RUNTIME sequence wait() 47: a QR-based solution which consists in estimating the norm of A -1 by computing the norm of R -1 with A = QR. This solution, which turns out to be less costly, alleviates the pivoting issue all together, uses only regular tile algorithms and allows code portability across various architectures, thanks to the underlying runtime system. The third section of the algorithm, rows 21 to 48, is the main loop of the algorithm, which iterates on U k and converges to the polar factors. This section of the algorithm is straightforward and follows the mathematical description of the problem using either a QR or a Cholesky factorization to calculate the next U . Finally, the last section, rows 49 to 53, computes the Hermitian polar factor H from the polar factor computed out of the main loop.

k = k +

C. Code Optimizations

The Chameleon library provides APIs to perform dense matrix computations. The first one, MORSE_xxxx_Tile, is a synchronous implementation of a linear algebra operation. This means that all the tasks required for the computations are submitted to the runtime, and then the library internally waits for the completion of all tasks before returning the control to the programmer. This is the first version we implemented, in black in the algorithm 1. To highlight the benefit of using task-based programming model (through tile algorithms) as opposed to the fork-join paradigm, as implemented in the LAPACK library, we have manually integrated synchronization points within the QR/Cholesky factorization kernel calls, at the end of each panel and update computations, to better emphasize on the performance discrepancy between both approaches. We refer to this reference implementation as Sync.

The second API, MORSE_xxxx_Tile_Async, ensures that all the tasks of an algorithm are submitted to the runtime, but their completion is not ensured when the function call returns. Thus, it is possible to simultaneously submit tasks of multiple operations. The runtime is in charge of keeping the data coherency of tasks, generated from different kernel calls, since these tasks may operate on the same data. Operations that are asynchronously submitted to the runtime are indicated in red in Algorithm 1. At some point of the algorithm, synchronization points are however required to guarantee the consistency of the results. This is made through a call to RUNTIME_sequence_wait(), which waits for the completion of all tasks. It is then possible to release synchronization in the three steps of the algorithm to ensure a better occupancy of the resources, especially on small to medium test cases, as presented in Section VIII. We refer to this implementation as Async. It is also noteworthy to mention that it is possible to estimate offline the minimal number of iterations performed in the main loop. In that case, the synchronization in line 45 can be safely removed for the first iterations and introduced only for the last iteration as a sanity check on the value conv against the convergence threshold.

The last optimization is the possible acceleration of the QRbased Halley iterations. This optimization consists in exploiting the identity matrix structure of the C 2 matrix in the QR factorization (line 34 in Algorithm 1) and the corresponding Q generation (line 35 in Algorithm 1). Indeed, thanks to tile algorithms, it is possible to design a specific QR factorization algorithm in order to factorize a dense matrix on top of an identity matrix. This new QR factorization takes into account the identity matrix structure so that only non-zeros tiles are operated on during the factorization. By the same token, during the Q generation step, only the non-zeros tiles containing the Householder reflectors will be accessed. This optimization is important as it reduces the number of flops as well as data movement. We refer to this implementation as OptId.

VII. ARITHMETIC COMPLEXITY

In this Section, we present the algorithmic complexity (flops) of the polar decomposition using two variants based on the Halley iteration (QDWH) and the SVD. For simplicity purposes, we consider only square dense matrices, but QDWH works also for rectangular matrices [START_REF] Nakatsukasa | Stable and Efficient Spectral Divide and Conquer Algorithms for the Symmetric Eigenvalue Decomposition and the SVD[END_REF].

A. The QDWH-based Polar Decomposition

The condition number estimation l 0 can be calculated using the LU factorization, which requires 2 3 n 3 , followed by two triangular solvers LX = Id and U A -1 = X, adding 2n 3 flops. Alternatively, l 0 can be calculated using the QR factorization, A = QR which needs 4 3 n 3 , followed by inverting the upper triangular matrix R with 1 3 n 3 . Calculating l 0 using the QR factorization needs less flops overall. Moreover, the resulting QR factors can be reused during the first iteration of QDWH, thanks to fine-grained computations.

As shown in Equation 1, the QDWH flops using QR-based iteration includes the QR decomposition of 2n × n matrix for a cost of (3 + 1

3 )n 3 flops. Then, forming

Q 1 Q 2 explicitly, needs (3 + 1 3
)n 3 flops. The product Q 1 Q 2 requires 2n 3 flops. Therefore, the arithmetic cost of each QR-based iteration is (8 + 23 )n 3 flops. For the Cholesky-based iteration in Equation 4, matrix-matrix multiplication involves 2n 3 , the Cholesky factorization needs 1 3 n 3 , and solving two linear systems requires 2n 3 . Therefore, the arithmetic cost of Cholesky-based iteration is (4 + 1

3 )n 3 . Computing the Hermitian polar factor H = U p A requires 2n 3 . Hence, the overall cost of QDWH is

(8+ 2 3 )n 3 ×#it QR +(4+ 1 3 )n 3 ×#it Chol +2n 3
, where #it QR and #it Chol correspond to the number of QR-based and Cholesky-based iterations, respectively. As discussed in [START_REF] Nakatsukasa | Stable and Efficient Spectral Divide and Conquer Algorithms for the Symmetric Eigenvalue Decomposition and the SVD[END_REF], the flop count of QDWH depends on l 0 , which involves during the QDWH iteration. The total flop count of QDWH for dense matrices ranges then from (10 + 2

3 )n 3 (for l 0 ≈ 1 with #it Chol = 2) to 41n 3 (for l 0 1, with typically #it QR = 3 and #it Chol = 3). Furthermore, taking advantage of the trailing identity matrix structure in the QR factorization (OptId) reduces the flop count of the iteration in Equation 1. Forming the upper triangular matrix R by applying the Householder reflectors with n + 1 nonzero elements Π n-1 k=1 H k A = R to k vectors requires 4(n + 1)k flops, therefore forming R needs Σ n-1 k=1 4(n -k)(n + 1) = 2n 3 flops. Accumulating the Householder reflectors to form

Q 1 Q 2 = Π n-1 k=1 H n-k I n 0 requires Σ n-1
k=1 4nk = 2n 3 , as explained in [START_REF] Nakatsukasa | Stable and Efficient Spectral Divide and Conquer Algorithms for the Symmetric Eigenvalue Decomposition and the SVD[END_REF]. Table I summarizes the total flop count of QDWH (including condition number estimation and Halley iteration) (1) when using LU to estimate l 0 (original implementation), (2) when using QR to estimate l 0 and reusing the QR factors in the first iteration of QDWH and (3) when additionally taking advantage of the identity matrix structure in QR-based iterations (Equation 1).

B. The SVD-based Polar Decomposition

The polar decomposition can be calculated via SVD as follows, A = U ΣV , then, U p = U V , H = V ΣV . Therefore, the flop count of this approach includes the cost of an SVD decomposition, a matrix-matrix multiplication to compute the orthogonal polar factor U p and a matrix-matrix well ill (1) QDWH+LU multiplication to calculate the Hermitian polar factor H. The standard approach to compute the SVD of a dense matrix is to first reduce it to bidiagonal form A = U 1 BV 1 . The subsequent left and right singular vectors from the bidiagonal solver are then accumulated during the back transformation phase, i.e., U = U 1 U 2 and V = V 2 V 1 , to calculate the singular vectors of the original matrix A. The final estimated flop count to calculate the SVD is 22n 3 , as implemented in the divide-and-conquer DGESDD [START_REF] Hansen | Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion[END_REF]. Then, we need to add 2n 3 to compute U p = U V , and n 3 to compute H = V ΣV (symmetric rank-k update operation). The final estimated cost of the polar decomposition using SVD is, therefore, 25n 3 . Compared to the QDWH-based polar decomposition (3) in Table I, this is 30% less than in case of ill-conditioned matrices and almost twice the flops in case of well-conditioned matrices. In theory, it seems there is a clear advantage to use SVD-based for the polar decomposition in presence of illconditioned matrices. However, the SVD algorithm inherently suffers from lack of parallelism, due to a very expensive panel factorization phase and may not be as competitive as QDWHbased approaches.

(13 + 1 )n 3 ≤ • • • ≤ (43 + 2 3 )n 3 (2) QDWH+QR (12 + 1 3 )n 3 ≤ • • • ≤ (41 + 1 3 )n 3 (3) QDWH+QR+OptId (12 + 1 3 )n 3 ≤ • • • ≤ (33 + 1 3 )n 3

VIII. PERFORMANCE RESULTS AND ANALYSIS

This Section provides a comprehensive performance analysis of the task-based QDWH implementation in the context of the Chameleon library with the dynamic runtime system StarPU on various architectures.

A. Environment Settings

We have considered three different systems, which are representative of the current manycore-based hardware trends. The first system is composed of dual-socket 16-core Intel Haswell Intel Xeon CPU E5-2698 v3 running at 2.30GHz equipped with 8 NVIDIA K80s (2 GPUs per board). The second system hosts the latest Intel commodity chip with dualsocket 14-core Intel Broadwell Intel Xeon E5-2680 v4 running at 2.4GHz. The third system has the latest Intel manycore Knights Landing (KNL) 7210 chips with 64 cores. For simplicity purposes, each system is named after his chip codename. Our QDWH implementation has been compiled with Intel compiler 16 and linked against the Chameleon library v0.9.0 with hwloc v1.11.4, StarPU v1.2.0 and Intel MKL v11.3.1. We have mainly considered ill-conditioned randomly generated matrices, since this represents the worse case scenario, where QDWH performs a maximum of six iterations. In particular, in the subsequent experiments, our QDWH implementation switches from Equation 1 to Equation 4 if c k is smaller than 100 (see Algorithm 1), which generates QR-based iterations for the first three followed by three Cholesky-based iterations.

B. Numerical Accuracy

We recall the polar decomposition of a given general matrix A ∈ R n×n : A = U p H. The norm . F denotes the Frobenius norm. To highlight the numerical robustness of the method, we use the following two accuracy metrics:

I-U p Up F A F
for the orthogonality of the polar factor U p and

A-UpH F A F
for the accuracy of the overall computed polar decomposition. Figure 1 presents the orthogonality of U p and the accuracy of the polar decomposition A = U p H for ill-conditioned matrix on the KNL system (very similar numerical results on other systems). We can distinguish two clusters, i.e., QDWHbased and SVD-base polar decomposition, with up to two digits difference in the orthogonality and accuracy magnitudes. Although both mostly employ orthogonal transformations, the SVD variant of the polar decomposition necessitates the QR algorithm, which may show some convergence limitations with ill-conditioned matrices, as shown later in Section VIII-F. Assessing the numerical accuracy/robustness of the task-based QDWH.

C. Incremental Optimizations

Figure 2 highlights the performance impact of various incremental optimizations on the task-based QDWH, as described in Section VI-C. Taking advantage of the identity matrix structure (OptId) engenders up to 20% performance improvements compared to the oblivious approach on all studied systems.

Running additionally in asynchronous mode (Async) further time to solution (up to 2.8x), especially for medium range of matrix sizes, where processing units run out of work and look-ahead techniques jump right in to fill the performance gap. For asymptotic matrix sizes, although work is abundant, the asynchronous mode still provides additional performance. In particular, on KNL and Haswell+8x80 systems, data movement engendered by NUMA and PCIe channels is expensive and can be overlapped by computations, thanks to the Async optimization.

D. Execution Traces

Figure 3 shows the execution traces when running in synchronous ( Tile API) and asynchronous ( Tile Async API) modes. We have added additional synchronization points within the Tile kernel API, after each panel/update computation, so that we can better capture the performance gain against coarse-grained computations engendered by block algorithms, as described in Section V. These traces have been obtained on the KNL system for a matrix size of 10K. Since the matrix is ill-conditioned, the task-based QDWH performs six iterations (three QR-based and three Cholesky-based). The green, blue and yellow blocks correspond to QR, Cholesky/Level 3 BLAS and Level 1/2 BLAS, respectively. We can clearly notice the idle time during the first three QR-based iterations when running with a synchronous mode (Figure 3 graphs, the performance curves of the task-based QDWH correspond to performance when all optimizations are enabled (i.e., Async and OptId).

E. Performance Scalability

Figure 4 demonstrates the performance scalability of the task-based QDWH implementation. The scalability is almost linear for the commodity CPU systems (i.e, Haswell/Broadwell). On KNL and Haswell+8xK80, although the overhead of moving data on these systems is higher than the commodity CPU platforms, StarPU is able to cope with these communication overheads and the overall scalability is still decent.

F. Performance Comparisons of QDWH Variants

Figure 5 reports task-based QDWH performance against other various existing QDWH implementations on ill (left) and well (right) conditioned matrices, across the three systems. Missing data points correspond to runs, which did not achieve the proper accuracy, as defined in Section VIII-B. For well-conditioned matrices, time to solution is much more shortened for the QDWH implementation variants, thanks to less iterations for convergence. The SVD variants of the polar decomposition do not seem to take advantage of such matrices since the bidiagonal reduction and the matrixmatrix multiplication have still to be performed in the same manner, regardless of the matrix condition number. All in all, the task-based QDWH achieves gains up to [6%, 8%] on Haswell and [39%, 17%] on Broadwell, [85%, 82%] on Haswell+8xK80, and [63%, 67%] on KNL against the best (non task-based) implementation for [ill, well]-conditioned matrices, respectively. Highest performance are achieved on systems where data movement are most expensive (e.g., NUMA for KNL and PCIe for Haswell+8xK80) since the asynchronous mode can mitigate the overhead of data transfers by overlapping communications with task computations. Also, compared to MAGMA QDWH [START_REF] Sukkari | A High Performance QDWH-SVD Solver Using Hardware Accelerators[END_REF], the task-based QDWH achieves gains up to [71%, 22%] on Haswell+4xK80 for [ill, well]-conditioned matrices, respectively.

G. Performance Comparisons Across Architectures

We have additionally considered two more recent architectures, i.e., a dual-socket 10-cores IBM Power8 (3.69GHz) and a dual-socket 16-cores Intel Haswell equipped with four NVIDIA Pascal P100 GPUs. Figure 6 presents the performance of the task-based QDWH across all systems investigated in the paper. The main idea is not to cross-compare the performance delivered by each system but rather to show that the task-based QDWH can support various architectures with a decent sustained peak (up to 90% and up to 60% of the sustained Chameleon DGEMM peak for CPU only systems and for KNL/GPUs platforms, respectively).

IX. CONCLUSION AND FUTURE WORK

We have presented a comprehensive performance analysis of a novel asynchronous task-based QDWH algorithm for the polar decomposition of a dense matrix. Thanks to finegrained computations, we have reduced by 20% the overall complexity by taking advantage of the identity structure of the matrix during the iterations, while exposing look-ahead opportunities to increase hardware occupancy. Furthermore, the Chameleon library and its dynamic runtime system StarPU abstracts the hardware complexity from end-users and is capable of asynchronously scheduling computational tasks on the underlying processing units. Thanks to its wide hardware range support, we demonstrated that StarPU can port a single sequential source code to a myriad of hardware systems. Experimental results of the asynchronous task-based QDWH show significant performance improvement (up to an order of magnitude) against state-of-the-art implementations on ill and well-conditioned matrices across various hardware technologies, which are paving the road to future petascale/exascale systems. Future works include using the task-based QDWH as a building block for the dense symmetric eigensolver and SVD on shared and distributed-memory systems. We would like also

Algorithm 1 1 5 2 6: 19 : 3 √

 112193 QDWH pseudo-code on top of Chameleon 1: /* Estimate the condition number */ 2: dlacpy Async( A, U ) U = A 3: dlacpy Async( A, B ) B = A 4: Anorm = dlange Async( A ) A : dgenm2( A, α ) α ≈ A RUNTIME sequence wait() 7: /* Compute U0 and l0*/ 8: dlascl Async( U, 1./α ) l0 = (α/1.1) * l0 20: 21: /* Compute the polar decomposition A = UpH using QDWH */ 22: k = 1, Li = β × α/1.1, conv = 100 23: while (conv ≥ 5eps || |Li -1| ≥ 5eps) do 24:

  Orthogonality of Up.

Fig. 1 .

 1 Fig. 1.Assessing the numerical accuracy/robustness of the task-based QDWH.

Fig. 2 .

 2 Fig.2. Assessing the performance of various incremental optimizations.

  Figure3shows the execution traces when running in synchronous ( Tile API) and asynchronous ( Tile Async API) modes. We have added additional synchronization points within the Tile kernel API, after each panel/update computation, so that we can better capture the performance gain against coarse-grained computations engendered by block algorithms, as described in Section V. These traces have been obtained on the KNL system for a matrix size of 10K. Since the matrix is ill-conditioned, the task-based QDWH performs six iterations (three QR-based and three Cholesky-based). The green, blue and yellow blocks correspond to QR, Cholesky/Level 3 BLAS and Level 1/2 BLAS, respectively. We can clearly notice the idle time during the first three QR-based iterations when running with a synchronous mode (Figure3(a)). The performance impact of synchronous execution for the next three Cholesky-based iterations is not as severe as QR-based iterations because the Cholesky panel factorization involves only the diagonal block (Figure3(b)). For the subsequent

Fig. 3 .

 3 Fig. 3. Assessing synchronous Vs asynchronous execution traces of taskbased QDWH on the KNL system with a matrix size of 10K.

Fig. 4 .

 4 Fig. 4. Assessing the task-based QDWH scalability.

ACKNOWLEDGMENT

The authors would like to thank Samuel Thibault from Inria for his support with StarPU, Jack Poulson from Google Inc. for his help in tuning Elemental and the vendors Cray/IBM/Intel/NVIDIA for their hardware donations and/or systems' remote accesses in the context of the Cray Center of Excellence, the Intel Parallel Computing Center and the NVIDIA GPU Research Center awarded to the Extreme Computing Research Center at KAUST.

 2 4 2 0 4 8 3 0 7 2 4 0 9 6 5 1 2 0 6 1 4 4 7 1 6 8 8 1 9 2 9 2 1 6 1 0 2 4 0 1 1 2 6 4 1 2 2 8 8 1 3 3 1 2 1 4 3 3 6 1 5 3 6 0 1 6 3 8 4 1 7 4 0 8 1 8 4 3 2 1 9 4 5 6 2 0 4 8 0 2 1 5 0 4 2 2 5 2 8 2 3 5 5 2 2 4 5 7 6 2 5 6 0 0 2 6 6 2to investigate other QDWH variants, which may require more flops but entails an even higher level of concurrency.


Mark Cummins

Michael Dowling

Fearghal Kearney

Oil Market Modelling: A Comparative Analysis of Fundamental and Latent Factor Approaches

Keywords: oil futures, fundamental models, latent factors, Vuong model comparison C32, G13, Q41

We formally compare fundamental factor and latent factor approaches to oil price modelling. Fundamental modelling has a long history in seeking to understand oil price movements, while latent factor modelling has a more recent and limited history, but has gained popularity in other financial markets. The two approaches, though competing, have not formally been compared as to effectiveness. For a range of short-medium-and long-dated WTI oil futures we test a recently proposed five-factor fundamental model and a Principal Component Analysis latent factor model. Our findings demonstrate that there is no discernible difference between the two techniques in a dynamic setting. We conclude that this infers some advantages in adopting the latent factor approach due to the difficulty in determining a well specified fundamental model.

Introduction

The modelling of prices in oil markets has received considerable attention in the literature. A particular focus of studies has been the identification of fundamental factors that drive oil markets. Everything from internal market characteristics, to linkages with other traded markets, to a wide range of macroeconomic factors, has been tested to determine an optimal fundamental model. In contrast to fundamental modelling a more recent development is based on identifying latent factors. These capture market dynamics using, for instance, principal component analysis of historical data in a purely statistical manner. Latent factors bypass the need to identify a fundamental model; especially useful given the debate as to precisely what fundamental variables should be in such a model. This paper sets out to determine whether fundamental modelling or the latent factor technique offers a better approach to modelling oil prices.

In terms of fundamental factors, as noted, many different lines of investigations have been pursued. Macroeconomic factors have been considered by [START_REF] Hamilton | Oil and the macroeconomy since world war ii[END_REF], [START_REF] Bernanke | Systematic monetary policy and the effects of oil price shocks[END_REF], [START_REF] Hamilton | What is an oil shock?[END_REF] and [START_REF] Barsky | Oil and the macroeconomy since the 1970s[END_REF]. How supply and demand imbalances, along with the movements in the business cycle, impact on oil prices has been investigated in depth [START_REF] Kilian | Not all oil price shocks are alike: Disentangling demand and supply shocks in the crude oil market[END_REF][START_REF] Kilian | The economic effects of energy price shocks[END_REF][START_REF] Askari | Oil price dynamics (2002-2006)[END_REF][START_REF] Kilian | Exogenous oil supply shocks: how big are they and how much do they matter for the us economy?[END_REF][START_REF] Kaufmann | The role of market fundamentals and speculation in recent price changes for crude oil[END_REF]. [START_REF] Heath | Macroeconomic factors in oil futures markets[END_REF] provides a contemporary study on the macroeconomic factors that drive oil prices, showing that measures of real economic activity forecast oil futures prices and, most notably, that real economy shocks have a resulting impact on oil prices.

The link between equity markets and oil markets has also been well studied, with the former representing a proxy for overall economic activity. [START_REF] Kilian | The impact of oil price shocks on the us stock market*[END_REF] show that US equity markets respond differently to demand and supply shocks in oil markets, and that in the aggregate both forms of shock account for a high percentage (22%) of the long-run variation in stock returns. [START_REF] Arouri | Does crude oil move stock markets in europe? a sector investigation[END_REF] examines the European stock markets and shows a similar link between the markets, although the strength of association is shown to vary significantly across sectors. [START_REF] Hammoudeh | Behavior of gcc stock markets and impacts of us oil and financial markets[END_REF] examine Gulf region stocks and show that, among other factors, oil price shocks influence equity prices. In the case of Oman and Saudi Arabia, oil prices are reported to account for up to 30% of the variation in stock returns. Other such studies that examine the interdependence between equity markets and oil markets include [START_REF] Hammoudeh | Dynamic relationships among gcc stock markets and nymex oil futures[END_REF], [START_REF] El-Sharif | Evidence on the nature and extent of the relationship between oil prices and equity values in the uk[END_REF], [START_REF] Filis | Dynamic correlation between stock market and oil prices: The case of oil-importing and oil-exporting countries[END_REF]. Several related studies investigate a similar question from the perspective of volatility transmission and spillover [START_REF] Malik | Shock and volatility transmission in the oil, us and gulf equity markets[END_REF][START_REF] Arouri | Does crude oil move stock markets in europe? a sector investigation[END_REF][START_REF] Creti | On the links between stock and commodity markets' volatility[END_REF].

Given that the global currency for oil trade is the USD, much research has focused on foreign exchange effects. [START_REF] Lizardo | Oil price fluctuations and us dollar exchange rates[END_REF] demonstrate that oil prices significantly explain movements in the value of the USD against major currencies, with particular evidence that rises in oil prices lead to an appreciation (depreciation) of net exporter (importer) country currencies against the USD. Reboredo (2012) finds the dependence between oil prices and exchange rates to be weak in general, although it is found that this dependence rose substantially in the aftermath of the global financial crisis.

Using a copula-GARCH approach, [START_REF] Aloui | Conditional dependence structure between oil prices and exchange rates: a copula-garch approach[END_REF] find rises in oil prices to be associated with the depreciation of the USD. Utilising alternative wavelet multi-resolution analysis, [START_REF] Reboredo | A wavelet decomposition approach to crude oil price and exchange rate dependence[END_REF] find that oil prices and exchange rates were not dependent in the pre-crisis period but that this changed to negative dependence with the onset of the crisis. Oil prices are shown to lead exchange rates and vice versa in the crisis period, but not in the pre-crisis period.

In addition to these identifiable fundamental factors, there is significant ongoing discourse over the increasing participation of financial players in the oil markets, and so the impact of speculation on both oil prices and volatility. The diversification benefits of commodities has been argued by many authors [START_REF] Bodie | Risk and return in commodity futures[END_REF][START_REF] Erb | The strategic and tactical value of commodity futures[END_REF][START_REF] Gorton | Facts and fantasies about commodity futures[END_REF] and as a result, the oil markets, and broader commodity markets, have seen increased financial player participation. The evidence on the effects of speculation is mixed. In some studies, speculation is shown to have a statistically significant effect on price and volatility in oil markets, particularly in the lead up to the historically high oil prices achieved pre-crisis [START_REF] Sornette | The 2006-2008 oil bubble: Evidence of speculation, and prediction[END_REF][START_REF] Kaufmann | Oil prices, speculation, and fundamentals: Interpreting causal relations among spot and futures prices[END_REF][START_REF] Cifarelli | Oil price dynamics and speculation: A multivariate financial approach[END_REF][START_REF] Du | Speculation and volatility spillover in the crude oil and agricultural commodity markets: A bayesian analysis[END_REF]. In contrast, [START_REF] Büyüksahin | Do speculators drive crude oil futures prices[END_REF] find no evidence that non-commercial positions, including hedge fund positions, have a causal effect on oil prices. [START_REF] Sanders | The adequacy of speculation in agricultural futures markets: Too much of a good thing?[END_REF], [START_REF] Irwin | Testing the masters hypothesis in commodity futures markets[END_REF], [START_REF] Irwin | Commodity index investment and food prices: does the masters hypothesis explain recent price spikes?[END_REF] similarly find no evidence for the influence of speculation across commodity markets, which is corroborated by the surveyed evidence reported by [START_REF] Irwin | Index funds, financialization, and commodity futures markets[END_REF].

In light of this wide and sometimes contradictory literature on the fundamental factors in oil markets, the real challenge with fundamental factor ap-proaches lies in the process of searching out and identifying appropriate variables and then testing their relevance in contemporaneously and dynamically modelling oil price movements. In contrast to fundamental factor approaches, latent factor approaches have been proposed as a purely statistical way of capturing market dynamics. Such methodologies circumvent the need for the explicit identification of fundamental factors. Latent factor approaches have been used to capture the complex movement in commodity forward curves [START_REF] Gibson | Stochastic convenience yield and the pricing of oil contingent claims[END_REF][START_REF] Schwartz | The stochastic behavior of commodity prices: Implications for valuation and hedging[END_REF][START_REF] Casassus | Stochastic convenience yield implied from commodity futures and interest rates[END_REF][START_REF] Trolle | Unspanned stochastic volatility and the pricing of commodity derivatives[END_REF][START_REF] Casassus | Economic linkages, relative scarcity, and commodity futures returns[END_REF][START_REF] Hamilton | Risk premia in crude oil futures prices[END_REF] and to forecast future oil prices [START_REF] Cabbibo | String theory[END_REF]Fiorenzani, 2004, Chantziara and[START_REF] Chantziara | Can the dynamics of the term structure of petroleum futures be forecasted? evidence from major markets[END_REF]; much of the the former literature focused on derivatives pricing and risk management applications. [START_REF] Chantziara | Can the dynamics of the term structure of petroleum futures be forecasted? evidence from major markets[END_REF] argue that principal component analysis (PCA) provides a parsimonious, non-parametric methodology to describe commodity forward curve dynamics, and should contain full information on the underlying economic drivers of oil prices. In contrast, [START_REF] Heath | Macroeconomic factors in oil futures markets[END_REF] argues that real economic activity has material effects on risk premiums and forecasts of oil futures prices, over and above the information in current futures prices.

Linking the two strands of modelling approach remains a gap in the literature to a large extent. [START_REF] Dempster | Determinants of oil futures prices and convenience yields[END_REF] consider the problem of linking latent factors to economic variables within the context of a proposed three-factor model for spot crude oil prices. Log spot prices are defined as an affine combination of latent short-, medium-and long-term factors, from which an analytic expression for log forward prices is derived as a function of all three factors. Using an EM algorithm combining Kalman filtering and MLE to extract out the latent factors, the authors formally test the relationships of these factors with a range of financial variables (including USD, S&P500, and VIX indices, as studied here), business cycle variables, demand variables, and trading variables. It is shown that the short-term factor links to the demand and trading variables, the medium-term factor links to the business cycle variables, and the long-term factor links to the financial variables. Thus there is some evidence of a link between latent factors and fundamental models as a technique to modelling, however as yet it has not been determined which approach is superior.

Motivated by the work of [START_REF] Dempster | Determinants of oil futures prices and convenience yields[END_REF], we make a number of contributions. In a first contribution, we conduct a formal model comparison of alternative fundamental factor and latent factor model specifications for oil futures returns. Specifically, for our fundamental factor model, we consider the five factors examined by [START_REF] Andreasson | Impact of speculation and economic uncertainty on commodity markets[END_REF]: the S&P 500 index; VIX volatility index; FRED trade-weighted USD currency index; US economic uncertainty index; and Working's T speculation index. For our latent factor model, we employ standard PCA methods to filter out the three most important principal components of the oil forward curve and use these as appropriate latent factors. Given the non-nested nature of the linear model specifications, we employ the non-nested model selection test of [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF] for our analysis. To the author's knowledge such a formal comparison of fundamental factor and latent factor models has not been conducted to date and so our study fills this gap.

In a second contribution, we conduct our model comparison exercise across the term structure of oil prices to examine whether there is a difference in model fit based on contract maturity. As noted already, [START_REF] Dempster | Determinants of oil futures prices and convenience yields[END_REF] argue a difference in the factors impacting on the short-, mediumand long-term. Our analysis across the forward curve provides insights into the ability of the respective fundamental and latent factor models to capture term structure movements.

Finally, in a third contribution, we explore the relevance of the speculation factor to the specification of the fundamental factor model for oil prices. For our analysis we consider the continuous monthly WTI futures CL1-CL9, CL12 and CL18 (corresponding to the 1-9, 12 and 18 month maturities respectively), with the sample spanning 10th January 2007 to 4th March 2016. Applying the non-nested model selection test of [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF] we find that for all contract maturities, the fundamental factor model, whether inclusive of speculation or not, and the latent factor model are equivalent in fit. This finding is relevant as latent factor based linear models have not received the same level of attention as fundamental factor based linear models for oil price modelling, yet there are some clear advantages: avoiding the need to search out and identify appropriate economic factors; avoiding criticisms about omitted variable and/or irrelevant variable bias; and being purely driven by data. As argued by [START_REF] Chantziara | Can the dynamics of the term structure of petroleum futures be forecasted? evidence from major markets[END_REF], the principal components should contain all of the information contained in the underlying economic variables that actually drive oil prices, and so should capture as much, if not more, information than a well specified fundamental factor model.

A secondary finding from our analysis is that speculation does not emerge as an important factor in the fundamental model specification, despite recent suggestions that it should be included in fundamental models. In our nonnested test involving the four-factor fundamental model, which excludes the speculation variable, we find that this reduced form model specification is also equivalent in terms of model fit to the principal components model. In our direct comparison of the five-factor and four-factor fundamental models, we find no evidence that the former model offers a superior model fit.

The remainder of the paper is organised as follows. Section 2 presents the fundamental factor and latent factor model specifications, along with the technical detail of the [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF] non-nested model selection test. Section 3 presents the results of our model comparison testing. Section 4 concludes.

Modelling and Testing Methodology

The fundamental factor model we consider is motivated by the analysis of [START_REF] Andreasson | Impact of speculation and economic uncertainty on commodity markets[END_REF]. The authors investigate economic, financial and speculation based variables for modelling a range of commodities, including crude oil. Excess speculation, USD exchange rates and S&P 500 equity index returns are all shown to have some form of causal relationship with the commodities considered, with crude oil showing additional relationships with the VIX volatility index and an economic uncertainty index. The economic uncertainty index is an innovative inclusion and is motivated by a set of studies that establish a link between oil price shocks and economic policy uncertainty (Kang andRatti, 2013a,b, Antonakakis et al., 2014). We use all five factors in our fundamental model specification.

As in [START_REF] Andreasson | Impact of speculation and economic uncertainty on commodity markets[END_REF], excess speculation is captured via the usual Working's T index (Working, 1960) defined as follows:

W T = 1 + SS HL + HS if HS ≥ HL W T = 1 + SL HL + HS if HS < HL
where SS (SL) is the number of short (long) positions held by speculators and HH (HL) is the number of short (long) positions held by hedgers.

The S&P 500 index is used here, as commonly done, as a barometer of economic activity in the US. The dynamic fundamental model for oil returns we consider is therefore specified as follows:

CL M t = β 0 + β SP 500 SP 500 t-1 + β V IX V IX t-1 + β U SD U SD t-1 +β EP U EP U t-1 + β Spec Spec t-1 + ε t ,
where CL denotes the log return of the continuous WTI crude oil (CL) contract of maturity M; SP 500 denotes the log return of the S&P 500 index;

V IX denotes the log change in the VIX volatility index; U SD denotes the log return of the trade weighted US dollar index; EP U denotes the log change in the economic policy uncertainty index for the US; and Spec is the change in the Working's T speculation index.

To examine whether speculation is indeed a driving factor of oil prices, we reduce the five-factor fundamental model down by means of excluding the speculation variable. We therefore specify the following nested four-factor fundamental model:

CL M t = β 0 + β SP 500 SP 500 t-1 + β V IX V IX t-1 + β U SD U SD t-1 +β EP U EP U t-1 + ε t .
Any deterioration in the fit of the four factor model relative to the five factor model, as measured against the latent factor model, would point to an important influence of speculation. While this provides indirect evidence on how well specified the two fundamental models are, we go further and conduct a direct model comparison of the models through the nested model version of the test of [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF].

In contrast to the above, latent factor approaches offer a purely statistical alternative to the objective of oil market modelling. A rationale for such methods is that they circumvent the need for the explicit identification of fundamental factors. PCA is one of the primary tools used for this purpose, providing abstract orthogonalised factors derived directly from the data and reproducing the original correlation matrix of the observed variables [START_REF] Chantziara | Can the dynamics of the term structure of petroleum futures be forecasted? evidence from major markets[END_REF]; in our case, the forward curve of oil prices. [START_REF] Chantziara | Can the dynamics of the term structure of petroleum futures be forecasted? evidence from major markets[END_REF] use PCA in this way and consider the potential of the extracted principal components for forecasting across four commodity markets, including the benchmark WTI and Brent crude oils.

Closely aligned to the [START_REF] Chantziara | Can the dynamics of the term structure of petroleum futures be forecasted? evidence from major markets[END_REF] method of PCA extraction, but given our focus on crude oil only, we propose the following dynamic latent factor model for oil returns:

CL M t = β 0 + β P C1 P C1 t-1 + β P C2 P C2 t-1 + β P C3 P C3 t-1 + ε t ,
where P C1, P C2, P C3 denote the first, second and third principal components extracted from the log returns of the futures contracts.

The next section describes the formal model comparison framework and the use of the Vuong test [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF].

Vuong Test of Non-Nested Models

Consider two non-nested (or indeed nested or overlapping) models F θ = {f (y | z; θ) ; θ ∈ Θ} and G γ = {g (y | z; γ) ; γ ∈ Γ}. [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF] tests the null hypothesis that the two models are equivalent, in the sense that

E 0 [log f (y | z; θ * )] = E 0 [log g (y | z; γ * )] ,
against the alternative hypothesis that F θ is better than G γ , in the sense that

E 0 [log f (y | z; θ * )] > E 0 [log g (y | z; γ * )] ,
or the alternative hypothesis that G γ is better than F θ , in the sense that

E 0 [log f (y | z; θ * )] < E 0 [log g (y | z; γ * )]
.

Here E 0 [•] is the expectations operator with respect to the true joint distribution of (y, z), and θ * and γ * are the pseudo-true values of θ and γ respectively.

The appropriate test statistic is the likelihood ratio (LR) for model F θ against model G γ defined as

LR n θn , γn ≡ L f n θn -L g n (γ n ) = n t=1 log f Y t | Z t ; θn g (Y t | Z t ; γn ) ,
where θn and γn are the maximum likelihood (ML) estimates of θ * and γ * respectively. The model selection test proposed by [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF] states that if F θ and G γ are strictly non-nested then

under H 0 :LR n θn , γn /ω n √ n -→ N (0, 1) under H f :LR n θn , γn /ω n √ n -→ +∞ under H g :LR n θn , γn /ω n √ n -→ -∞
where H 0 is the null hypothesis that F θ and G γ are equivalent

H 0 : E 0 log f (Y t | Z t ; θ * ) g (Y t | Z t ; γ * ) = 0 tested against the alternative hypotheses H f that F θ is better than G γ H f : E 0 log f (Y t | Z t ; θ * ) g (Y t | Z t ; γ * ) > 0 or H g that F θ is worse than G γ H g : E 0 log f (Y t | Z t ; θ * ) g (Y t | Z t ; γ * ) < 0.
ω2 n is the sample analog of the log-likelihood variance ω 2 * ≡ var 0 log f (Yt|Zt;θ * ) g(Yt|Zt;γ * ) defined as follows:

ω2 n ≡ 1 n n t=1   log f Y t | Z t ; θn g (Y t | Z t ; γn )   2 -   1 n n t=1 log f Y t | Z t ; θn g (Y t | Z t ; γn )   2 ,
which is used in the variance test of [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF] to discriminate between two models if not equivalent in the sense that f

(• | •; θ * ) = g (• | •; γ * ). The associated hypotheses are H ω 0 : ω 2 * = 0 versus H ω A : ω 2 * = 0.
The implementation of the [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF] test is quite straightforward involving selecting an appropriate critical value c from the standard normal distribution. When the test statistic LR n θn , γn /ω n √ n is higher than c then F θ is deemed the better model, when it is lower than -c then G γ is deemed the better model, and when it lies in the intermediate region one cannot discriminate between the two models.

Empirical Results

For our analysis, we use the continuous monthly contracts CL1-CL9, CL12 and CL18 obtained from the CME. These contracts span the 1-9, 12 and 18 month maturities and allow us to consider a range of short-, mediumand long-term movements in the WTI forward curve. For the model com-parison exercise to follow, this gives sufficient flexibility to examine whether there is a difference in model fit based on contract maturity. The sample covers the period 10th January 2007 to 4th March 2016, observed daily. Table 1 provides summary statistics for the CL log returns series. Table 2 provides summary statistics for the five fundamental variables considered in our study. SP 500 is the S&P 500 index, V IX is the VIX volatility index, U SD is the trade weighted US dollar index, EP U is the economic policy uncertianty index for the US, and Spec is the speculation variable, defined here as the Working's T index (see Section 2). ∆ denotes daily log return/log change.

Table 3 and 4 present the regression results for the five-factor and four-factor fundamental models respectively. Notable in the former case is that speculation is not found to be statistically significant for any of the maturities considered. Economic uncertainty likewise is also found not to have a statistically significant effect on next period oil returns across the term structure.

When speculation is dropped in the four-factor case, no effect from economic uncertainty is again observed. The speculation finding aligns with the linear causality testing of [START_REF] Andreasson | Impact of speculation and economic uncertainty on commodity markets[END_REF], who find no causality in either direction between speculation and oil returns. In contrast though to our finding on economic uncertainty, [START_REF] Andreasson | Impact of speculation and economic uncertainty on commodity markets[END_REF] confirm a causal effect of this variable on oil returns. However, our result is found to be consistent for all maturities considered.

Broad equity performance, as a proxy for overall economic activity, can be seen to be a consistent driving factor of oil prices across the contract maturities, having a positive effect. Market volatility and USD currency exposure are also found to be statistically significant, although the positive effects reported are not observed across the entire forward curve. In fact, up to eight months both factors play a role, while the effects thereafter are only seen for the 12-month contract maturity. For the medium-term CL9 and the long-term CL18 contracts neither factor seems to emerge as having an influence.

For the latent factor model, Table 5 first provides the variance statistics for the first three principal components extracted from the forward returns data. The first principal component explains 96.4% of the variance, with the cumulative proportion rising to 99% with the inclusion of the second principal component and to 99.6% with the further inclusion of the third principal component. Table 6 reports the regressions results for the latent factor model specification based off these three estimated principal components.

Notable in the findings is that the first principal component emerges as being a statistically significant driver of oil returns for all maturities. Only for CL2 does the second principal component appear to be relevant, while only for CL1 does the third principal component display significance. This contrasts starkly with the analysis of [START_REF] Chantziara | Can the dynamics of the term structure of petroleum futures be forecasted? evidence from major markets[END_REF].

While their analysis is broader in scope, in that a more expansive twelvelatent factor model is considered using the first three principal components extracted from each of WTI, Brent, (NYMEX) heating oil and (NYMEX) gasoline futures returns, the authors report that it is only the third principal component derived from Brent returns, rather than WTI returns, that explains movements in WTI, and this is only evidenced for the intermediate CL3-CL5 contracts. The effect of the first principal component in our case is consistently observed across all maturities.

We proceed now to the formal model comparison exercise using the nonnested test of [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF]. 3 For completeness, we initially implement the variance test as proposed by [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF] to confirm that the two model specifications are not equivalent in the conditional density sense as outlined in Section 2.1. The results of the variance and non-nested model selection tests are presented in Table 7.

From the variance test results it can be seen that, whether comparing the principal component model against the five-factor or four-factor fundamental model, the competing models are deemed to be distinguishable, in the sense that the respective conditional densities are not equivalent; this is found at the 5% or 10% significance level across contract maturity, although in the latter case the p-values are marginally outside the 5% threshold. This justifies and prompts our use of the [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF] non-nested test of model selection in the next step. In all cases, without exception, the latent factor model and fundamental factor model, whether inclusive of speculation or not, appear to have equal fit to the population. The reported p-values exceed the 10% threshold and so we fail to reject the null hypothesis that the models are equivalent.

To explore the robustness of this finding, and in an attempt to discriminate between the two models, we generate confidence intervals for the model AIC and BIC statistics that are consistent with the [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF] theory.

The 95% confidence intervals are presented in Tables 8 and9; the former involving the five-factor fundamental model, the latter involving the fourfactor fundamental model. Here again it can be observed that the models are not distinguishable in terms of fit, and so one cannot be selected over the other in each pairwise case. We therefore conclude with confidence that there is no discernible distinction between our fundamental factor and latent factor specification in capturing movement in WTI crude oil prices.

This model comparison finding is significant given the fundamentally different approaches taken to modelling oil prices. In particular, while latent factor based linear models have not received the same level of attention as fundamental factor based linear models for oil price modelling, there as some clear advantages: avoiding the need to search out and identify appropriate economic factors; in so doing, avoiding criticisms about omitted variable and/or irrelevant variable bias; and being purely driven by data. In principle, the principal components should contain all of the information contained in the underlying economic variables [START_REF] Chantziara | Can the dynamics of the term structure of petroleum futures be forecasted? evidence from major markets[END_REF] that drive oil prices, and so should in theory capture as much, if not more, information than a well specified fundamental factor model.

As a secondary finding, it is notable that there is no notable distinction between the five-factor fundamental model and the four-factor counterpart, which excludes the speculation variable. This gives some indirect evidence to the absence of a role for speculation in driving oil prices. Given that this evidence is indirect in nature, we complete our analysis in this section with a direct application of the nested [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF] test to the two model specifications. Table 10 presents the results. The variance test shows that the five-factor and four-factor models are distinguishable, in the sense that the models respective conditional densities are not equivalent. However, the nested test results show that, for almost all maturities, we fail to reject the null hypothesis that the five-factor and four-factor models are indistinguishable. Only in the case of the 18-month maturity contract is the five-factor model deemed a better fit; this being concluded at the 10% significance level.

So we conclude overall that speculation does not have a significant influence on oil futures prices. Regression results for the five-factor fundamental model specified as follows:

CL M t = β 0 + β SP 500 SP 500 t-1 + β V IX V IX t-1 + β U SD U SD t-1 + β EP U EP U t-1 +
β Spec Spec t-1 + ε t where CL denotes the log return of the continuous WTI crude oil (CL) contract of maturity M; SP 500 denotes the log return of the S&P 500 index;

V IX denotes the log change in the VIX volatility index; U SD denotes the log return of the trade weighted US dollar index; EP U denotes the log change in the economic policy uncertainty index for the US; and Spec is the change in the speculation variable, defined here as the Working's T index (see Section 2). Standard errors in parenthesis. * * * , * * , * denote statistical significance at 1%, 5% and 10% respectively.

M t = β 0 + β SP 500 SP 500 t-1 + β V IX V IX t-1 + β U SD U SD t-1 + β EP U EP U t-1 + ε t where
CL denotes the log return of the continuous WTI crude oil (CL) contract of maturity M; SP 500 denotes the log return of the S&P 500 index; V IX denotes the log change in the VIX volatility index; U SD denotes the log return of the trade weighted US dollar index; and EP U denotes the log change in the economic policy uncertainty index for the US (see Section 2). Standard errors in parenthesis. * * * , * * , * denote statistical significance at 1%, 5% and 10% respectively. [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF]. The variance test is used to discriminate between two models, F θ and G γ , if not equivalent in the sense that f (

• | •; θ * ) = g (• | •; γ * ). The associated hypotheses are H ω 0 : ω 2 * = 0 versus H ω A : ω 2 * = 0
, where ω * is the log-likelihood variance. See Section 2.1 for details. The non-nested model selection test tests the null hypothesis that F θ and G γ are equivalent, H 0 : E 0 log f (Yt|Zt;θ * ) g(Yt|Zt;γ * ) = 0, against the alternative hypotheses that F θ is better than G γ , H f : E 0 log f (Yt|Zt;θ * ) g(Yt|Zt;γ * ) > 0, or that F θ is worse than G γ , H g : E 0 log f (Yt|Zt;θ * ) g(Yt|Zt;γ * ) < 0. See Section 2.1 for details. We designate F θ to be the latent (PC) factor model and G γ to be either the five-factor (5F) or four-factor (4F) fundamental model. P-values in parenthesis. that G γ fits as well as F θ against the alternative hypothesis H nested A that F θ fits better than G γ . We designate F θ to be the five-factor (5F) fundamental model and G γ to be the four-factor (4F) fundamental model. P-values in parenthesis.

Conclusion

In this study, we conduct a formal model comparison of fundamental factor and latent factor modelling approaches to capturing oil market movements. To the author's knowledge such a formal comparison has not been conducted to date and so our study fills this gap. We additionally contribute by means of conducting our model comparison exercise across the term structure of oil prices to examine whether there is a difference in model fit based on contract maturity. Prior literature [START_REF] Dempster | Determinants of oil futures prices and convenience yields[END_REF] suggests that different factors impact on the short-, medium-and long-term, and so our analysis across the forward curve provides insights into the ability of the respective fundamental and latent factor models to capture these term structure movements. In a final contribution, we use our model comparison setting to formally test whether speculation is a fundamental factor affecting oil prices.

We find the fundamental factor model and latent factor model specifications to be of equal fit to the overall population of oil futures prices. This finding is significant given the fundamentally different approaches taken to modelling oil prices. Latent factor models, which have received much less attention in the literature, outside of derivatives applications, offer real merit it seems relative to the popular approach of fundamental modelling. In particular, latent factor models avoid the need to search out and identify appropriate economic factors, in so doing avoid criticisms about omitted variable and/or irrelevant variable bias, and are accessible and easily implemented being purely driven by the data of interest. In principle, the principal com-ponents should contain all of the information contained in a well specified fundamental factor model.

The novelty of our study is the formal model comparison of fundamental factor and latent factor approaches. The work should motivate further research to appraise the two alternative strands of modelling, using an extended range of fundamental and latent factors. This should be performed on a more expansive cross-commodity basis to see if the findings observed here hold across markets.

  Our model comparison set up facilitates an examination of speculation in two ways. Firstly, we drop the speculation variable from the five factor fundamental model and in this four factor setting we re-test against the principal components based model. Any deterioration in the fundamental model fit relative to the principal components model would point to the influence of speculation. Secondly, and more formally, we use the nested model version of the Vuong (1989) test to perform a direct model comparison of the five factor and four factor fundamental models.

Table 1 :

 1 Descriptive Statistics: WTI CL Log-Return Series

	Mean	Std Dev Skew	Kurt
	∆CL1 0.0043%	0.0025	0.3022 5.6212
	∆CL2 0.0025%	0.0024	0.0356 3.7477
	∆CL3 0.0016%	0.0022 -0.0085 3.4223
	∆CL4 0.0013%	0.0022 -0.0113 3.2885
	∆CL5 0.0010%	0.0021 -0.0362 3.2501
	∆CL6 0.0019%	0.0021 -0.0232 3.2661
	∆CL7 0.0014%	0.0020 -0.0172 3.2785
	∆CL8 0.0009%	0.0020 -0.0207 3.2780
	∆CL9 0.0006%	0.0019 -0.0181 3.3825
	∆CL12 0.0014%	0.0018 -0.0301 3.3491
	∆CL18 -0.0011% 0.0017 -0.0736 3.6792
	CL1-CL9, CL12 and CL18 are the WTI continuous monthly futures contracts with corre-
	sponding 1-9, 12 and 18 month maturities respectively. ∆ denotes daily log return.

Table 2 :

 2 Descriptive Statistics: Fundamental Variables

		Mean	Std Dev Skew	Kurt
	∆SP 500 0.00019% 0.0134 -0.3554 8.0117
	∆V IX 0.00309% 0.0773	1.3242 5.9370
	∆U SD 0.00006% 0.0049 -0.3509 4.3325
	∆EP U 0.06711% 0.5855	3.3724 28.2354
	Spec	1.11873% 0.0361	0.7540 0.7308

Table 3 :

 3 Five-Factor Fundamental Model: Regression Results

	β 0	β SP 500	β V IX	β U SD	β EP U	β Spec
	CL1 -0.0177 0.2504 * * * 0.0202 * 0.4125 * * * 0.0005	0.0157
	(0.0165) (0.0610) (0.0104) (0.1099) (0.0009) (0.0147)
	CL2 -0.0149 0.2031 * * * 0.0165 * 0.2953 * * * -0.0001	0.0132
	(0.0152) (0.0564) (0.0096) (0.1016) (0.0008) (0.0136)
	CL3 -0.0139 0.2040 * * * 0.0152	0.2484 * * -0.0002	0.0124
	(0.0146) (0.0539) (0.0092) (0.0972) (0.0008) (0.0130)
	CL4 -0.0129 0.2048 * * * 0.0151 * 0.2105 * * -0.0003	0.0115
	(0.0146) (0.0539) (0.0092) (0.0972) (0.0008) (0.0130)
	CL5 -0.0116 0.2072 * * * 0.0145 * 0.1953 * * -0.0003	0.0104
	(0.0136) (0.0504) (0.0086) (0.0909) (0.0008) (0.0122)
	CL6 -0.0106 0.2060 * * * 0.0143 *	0.1716 *	-0.0004	0.0094
	(0.0133) (0.0491) (0.0084) (0.0884) (0.0007) (0.0118)
	CL7 -0.0097 0.2075 * * * 0.0139 * 0.1704 * * -0.0004	0.0084
	(0.0130) (0.0480) (0.0082) (0.0865) (0.0007) (0.0116)
	CL8 -0.0091 0.2124 * * * 0.0160 * * 0.1517 *	-0.0003	0.0081
	(0.0127) (0.0470) (0.0080) (0.0847) (0.0007) (0.0113)
	CL9 -0.0083 0.1725 * * * 0.0107	0.1359	-0.0002	0.0074
	(0.0124) (0.0459) (0.0079) (0.0828) (0.0007) (0.0111)
	CL12 -0.0065 0.2032 * * * 0.0138 *	0.1357 *	-0.0002	0.0058
	(0.0117) (0.0433) (0.0074) (0.0780) (0.0006) (0.0104)
	CL18 -0.0032 0.1897 * * * 0.0107	0.1157	-0.0003	0.0028
	(0.0107) (0.0397) (0.0068) (0.0715) (0.0006) (0.0096)

Table 4 :

 4 Four-Factor Fundamental Model: Regression Results Regression results for the four-factor fundamental model specified as follows: CL

	β 0	β SP 500	β V IX	β U SD	β EP U
	CL1 -0.0001 0.2482 * * * 0.0198 * 0.4117 * * * 0.0004
	(0.0005) (0.0609) (0.0104) (0.1099) (0.0009)
	CL2 -0.0001 0.2013 * * * 0.0162 * 0.2946 * * * -0.0001
	(0.0005) (0.0563) (0.0096) (0.1016) (0.0008)
	CL3 -0.0001 0.2023 * * * 0.0149	0.2477 * * -0.0002
	(0.0005) (0.0539) (0.0092) (0.0972) (0.0008)
	CL4 -0.0001 0.2032 * * * 0.0149 * 0.2099 * * -0.0003
	(0.0005) (0.0521) (0.0089) (0.0940) (0.0008)
	CL5 -0.0001 0.2057 * * * 0.0143 * 0.1948 * * -0.0003
	(0.0004) (0.0504) (0.0086) (0.0909) (0.0007)
	CL6 -0.0000 0.2047 * * * 0.0141 *	0.1711 *	-0.0005
	(0.0004) (0.0490) (0.0084) (0.0884) (0.0007)
	CL7 -0.0000 0.2063 * * * 0.0137 * 0.1700 * * -0.0005
	(0.0004) (0.0480) (0.0082) (0.0865) (0.0007)
	CL8 -0.0001 0.2113 * * * 0.0158 * * 0.1512 *	-0.0003
	(0.0004) (0.0470) (0.0080) (0.0847) (0.0007)
	CL9 -0.0001 0.1715 * * * 0.0105	0.1355	-0.0002
	(0.0004) (0.0459) (0.0079) (0.0828) (0.0007)
	CL12 -0.0001 0.2024 * * * 0.0136 *	0.1354 *	-0.0003
	(0.0004) (0.0432) (0.0074) (0.0780) (0.0006)
	CL18 -0.0001 0.1893 * * * 0.0107	0.1155	-0.0003
	(0.0003) (0.0396) (0.0068) (0.0715) (0.0006)

Table 7 :

 7 Vuong Test ResultsResults for the variance test and non-nested model selection test of

		F θ :PC v G θ :5F	F θ :PC v G θ :4F
		H ω 0 ;H ω A	H 0 ;H f H 0 ;H g H ω 0 ;H ω A	H 0 ;H f H 0 ;H g
	CL1	0.04	-0.46	-0.07	-0.41
		(0.051) (0.676) (0.324) (0.050) (0.658) (0.342)
	CL2	-0.02	-0.09	0.24	-0.04
		(0.047) (0.534) (0.466) (0.046) (0.516) (0.484)
	CL3	-0.48	-0.28	-0.36	-0.24
		(0.044) (0.610) (0.390) (0.044) (0.594) (0.406)
	CL4	-0.21	-0.36	-0.55	-0.32
		(0.045) (0.640) (0.360) (0.044) (0.625) (0.375)
	CL5	0.05	-0.50	0.31	-0.47
		(0.048) (0.691) (0.309) (0.047) (0.680) (0.320)
	CL6	-0.07	-0.61	0.18	-0.58
		(0.047) (0.729) (0.271) (0.047) (0.720) (0.280)
	CL7	-0.24	-0.61	-0.47	-0.59
		(0.052) (0.730) (0.270) (0.051) (0.723) (0.277)
	CL8	-0.07	-0.65	0.19	-0.63
		(0.053) (0.743) (0.257) (0.052) (0.737) (0.263)
	CL9	0.26	-0.45	-0.44	-0.43
		(0.041) (0.675) (0.325) (0.041) (0.668) (0.332)
	CL12	0.23	-0.87	0.33	-0.86
		(0.056) (0.807) (0.193) (0.056) (0.804) (0.196)
	CL18	0.36	-1.10	0.19	-1.10
		(0.058) (0.864) (0.136) (0.058) (0.864) (0.136)

Table 10 :

 10 Vuong Nested Test Results: Five-Factor v Four-Factor Fundamental Models Results for the variance test and nested model selection test of[START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF]. The variance test is used to discriminate between two models, F θ and G γ , if not equivalent in the sense that f (• | •; θ * ) = g (• | •; γ * ). The associated hypotheses are H ω 0 : ω 2 * = 0 versus H ω A : ω 2 * = 0,where ω * is the log-likelihood variance. See Section 2.1 for details. The nested model selection test tests the null hypothesis H nested

		F θ :5F v G γ :4F
		H ω 0 ;H ω A	H nested 0	;H nested A
	CL1	-0.07	0.23
		(0.00)	(1.00)
	CL2	-0.18	-0.58
		(0.00)	(0.95)
	CL3	0.72	-0.07
		(0.00)	(0.91)
	CL4	0.13	0.01
		(0.00)	(0.83)
	CL5	0.01	-0.28
		(0.00)	(0.72)
	CL6	-0.12	0.02
		(0.00)	(0.63)
	CL7	-0.12	0.10
		(0.00)	(0.56)
	CL8	-0.34	0.00
		(0.00)	(0.51)
	CL9	-0.01	-0.05
		(0.00)	(0.44)
	CL12	0.60	0.25
		(0.00)	(0.30)
	CL18 -0.07	0.00
		(0.00)	(0.08)
				0

The 10 newspaper sources are USA Today, Miami Herald, Chicago Tribune, Washington Post, Los Angeles Times, Boston Globe, San Francisco Chronicle, Dallas Morning News, New York Times, and Wall Street Journal.

2 www.policyuncertainty.com

The[START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF] test is implemented in R using the nonnest2 package. The package guide is available at https://cran.r-project.org/web/packages/nonnest2/nonnest2.pdf. Supporting documentation describing implementation and interpretation is available at https://cran.r-project.org/web/packages/nonnest2/vignettes/nonnest2.pdf.

0.0000 -0.0177 * * * 0.0233 0.0949 (0.0004) (0.0059) (0.0354) (0.0771) CL12 0.0000 -0.0158 * * * 0.0085 0.0732 (0.0003) (0.0056) (0.0334) (0.0728) CL18 -0.0000 -0.0132 * * * -0.0081 0.0199 (0.0003) (0.0051) (0.0306) (0.0668)

Regression results for the latent factor model specified as follows:

CL M t = β 0 + β P C1 P C1 t-1 + β P C2 P C2 t-1 + β P C3 P C3 t-1 + ε t , where P C1, P C2, P C3 denote the first, second and third principal components extracted from the log returns of the WTI futures contracts (see Section 2). Standard errors in parenthesis. * * * , * * , * denote statistical significance at 1%, 5% and 10% respectively. Working, H. (1960). Speculation on hedging markets. Food Research Institute Studies, (02).


Didier Hantz

Thomas Dewez

Clara Levy

Antoine Guerin

Michel Jaboyedoff

HANTZ Didier -Rockfall frequency in different geomorphological conditions Rockfall frequency in different geomorphological conditions

Keywords: rockfall, frequency, causal factors

The knowledge of the rockfall frequency in a cliff is needed for quantitative rockfall hazard assessment [START_REF] Hantz | A new approach of diffuse rockfall hazard[END_REF]. It can be estimated from historical data bases (e.g. [START_REF] Hantz | Rock fall hazard assessment: from qualitative to quantitative failure probability[END_REF] or from diachronic comparison of digital cliff topographic models (e.g. [START_REF] Dewez | Probabilistic coastal cliff collapse hazard from repeated terrestrial laser surveys : case study from Mesnil Val (Normandy, northern France)[END_REF][START_REF] Guerin | Investigating rock fall frequency using Terrestrial Laser Scanner[END_REF]). An empirical approach is proposed to be used when historical data bases are not significant (not enough events) and diachronic digital models not available.

VOLUME-FREQUENCY RELATION AND EROSION RATE

It is well recognised that the frequency of the rockfalls occurring in a rock wall decreases when the volume considered increases. It has appeared that the best law for describing this decrease is a power law, which can be applied either to the temporal frequency or to the spatial-temporal frequency (number of rockfalls per year and per m2 of cliff):

F = A V -B ( 1
)
Where F is the frequency of rockfalls bigger than V, A is the frequency of rockfalls bigger than 1 m3 , and B is the scaling exponent, whose value is usually between 0.4 and 0.8. Note that the relation holds up to a finite value of V because there is a maximal rockfall volume (V max ) which depends on the size of the cliff. This maximal volume is rarely observed because its return period is usually larger than the observation period. Integrating the volume V from 0 to V max (assuming B < 1) gives the volumetric retreat rate of the cliff [START_REF] Hantz | Rock fall hazard assessment: from qualitative to quantitative failure probability[END_REF]:

W = V max (1-B) A / (1-B) (2)
If the spatio-temporal frequency is used (F st = A st V -B ), Equation ( 2) gives the linear retreat rate (or erosion rate) of the cliff (E = V max (1-B) A st / (1-B)). As this rate includes the biggest possible volumes, it is usually higher than the rate obtained by summing the observed volumes or measuring the retreat of the cliff crest. V max is difficult to estimate, but one can note that if its uncertainty is a factor 10, the uncertainty on the erosion rate is lower (factor 10 (1-B) ). The retreat rate given by Equation (2) includes very large rockfalls (flow like movements) as well as smaller rockfalls (with negligible interaction between individual particles). As the methods used to simulate the propagation of these phenomena are different, it is relevant to calculate different volumetric retreat rates using this relation [START_REF] Hantz | Rock fall hazard assessment: from qualitative to quantitative failure probability[END_REF]:

W = [V 2 (1-B) -V 1 (1-B) ] AB / (1-B) ( 3
)
Where W is the volumetric retreat rate due to rockfalls whose volume is between V 1 and V 2 . A method to estimate the frequency of impacts due to small rockfalls is given by [START_REF] Hantz | A new approach of diffuse rockfall hazard[END_REF].

ESTIMATING THE ROCKFALL FREQUENCY PARAMETERS

Three parameters are necessary to characterize the rockfall activity of a cliff: V max , B and E or A st . V max and B are geometrical parameters. It is reasonable to assume they depend essentially of the rock mass structure (notably discontinuities extension and spacing). A st (or E) is a temporal parameter (year -1 ) which reflects the activity of the processes leading to rockfalls. The processes leading to landslides (of all types) have been listed by [START_REF] Popescu | A suggested method for reporting landslide causes[END_REF].

They can be divided in 3 groups: geomorphological processes, physical and chemical processes and human processes. [START_REF] Effendiantz | Les études spécifiques d'aléa lié aux éboulements rocheux[END_REF] have listed the processes leading to failure in rock walls. The more relevant ones have to be considered for estimating the erosion rate. Note that the Slope Mass Rating system [START_REF] Romana | New adjustment ratings for application of Bieniawski classification to slopes[END_REF] doesn't seem adequate for estimating the erosion rate, because it doesn't consider most of these processes (except groundwater action). It is useful to predict the stability of a rock cut rather than the evolution of an existing rock slope.

The rockfall frequency parameters B, A st and E have been determined for different types of rock masses, geomorphological and climatic conditions (Table 1). It seems that B is lower for massive rocks (massive limestone, chalk, gneiss) than for bedded rocks (bedded limestone). This finding is in agreement with the results obtained by [START_REF] Hungr | Magnitude and frequency of rock falls and rock slides along the main transportation corridors of southwestern British Columbia[END_REF] who concluded that a B value in the order of 0.7 is characteristic in moderately to highly jointed metamorphic, igneous and strong sedimentary rocks and that a lower value, of the order of 0.4, appears appropriate for massive felsic intrusive rocks (which possibly produce a relatively greater proportion of large-magnitude, structurally controlled failures). The erosion rate (and A st ) is strongly influenced by the geomorphological context: Erosion rates of some cm/year to some dm/year are obtained for coastal cliffs submitted to the waves action, while values of some tenth of mm/year to some cm/year are derived for middle mountain cliffs.

CONCLUSION

Orders of magnitude have been proposed for the rockfall frequency parameters according to the rock mass properties and geomorphological context.

  3rd RSS Rock Slope Stability conference, Lyon 2016 HANTZ Didier -Rockfall frequency in different geomorphological conditions

Table 1 : Rockfall frequency parameters

 1

	Rock	Massive	Massive	Bedded	Bedded	Massive	Massive
		limestone	limestone	limestone	limestone	chalk	gneiss
		(Urgonian)	(Urgonian)	(Valaginian)	(Sequanian)
	Site	Isère	Haute-Savoie	Gorgette	Saint-Eynard	Mesnil Val	Venosc
	Cliff area (m 2 )	6.1 10 7	7.8 10 7	5.1 10 4	1.3 10 5	3.7 10 4	3.5 10 5
	Period length (year)	62	22	3.2	3.2	2.3	3.2
	Number of rockfalls	87	45	147	344	8567	12
	B	0.52	0.60	0.57	0.75	0.54	0.38
	Ast (m -2 .year -1 )	1.4 10 -7	2.2 10 -7	9.6 10 -5	9.8 10 -5	1.1 10 -3	8.4 10 -6
	Maximal volume (m 3 )	10 7	10 7	10 6	10 6	10 5	10 6
	Erosion rate (m -1 .year -1 )	0.0007	0.0004	0.085	0.012	0.48	0.07

D'AMATO Julie, HANTZ Didier, ROSSETTI Jean-Pierre, BAILLET Laurent, ISTerre, Univ. Grenoble-Alpes, F, didier.hantz@ujf-grenoble.fr

DEWEZ Thomas, LEVY Clara, BRGM, F

GUERIN Antoine, JABOYEDOFF Michel, CRET, Faculty of Geosciences and Environment, University of Lausanne, CH

Acknowledgements. The authors thank the national project C2ROP and the Federation VOR for fundings.