Thomas Vallée
email: thomas.vallee@univ-nantes.fr

Computing the zeros and weights of Gauss-Laguerre and Gauss-Hermite quadratures: 2 Matlab files

This article provides two simple Matlab files to compute the zeros and weights of the Gauss-Laguerre and Gauss-Hermite quadratures.

Introduction

Economic problems usually deal with integrals of the form

∞ -∞ e -x 2 f (x)dx or ∞ 0 e -x f (x)dx.
The first one arises naturally in problem evolving normal random variables. The second one is simply a typical exponentially discounted sum. To evaluate (numerically) such integrals, one has to use appropriated methods, that is the Gauss-Hermite and the Gauss-Laguerre quadratures respectively.

These both procedures require to redefine these integrals as some special polynomials that are evaluated as the sum depending on some weights w i and zeros x i given the number of points to be used n. Although tables of these values exist [START_REF] Stroud | Gaussian Quadrature Formulas[END_REF], we provide two simple Matlab programs that calculate, given n, the optimal weights and zeros. Other methods to compute such weights and zeros [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF] exist and are based on some recursive procedure. The method we used is based on the fact that the associated polynomial can be rewritten as some finite series. and where w i are some weights defined by

w i = (n!) 2 x i (L n (x i )) 2 . 1
After some algebras, one can find that the polynomial L n (x) can be redefined as follows :

L n (x) = e x d n dx n (e -x x n ) (2.1) = n+1 i=1 (-1) n+1-i [b i n! (n + 1 -i)! x n+1-i ]
with 0! ≡ 1 and where the b i are some binomial coefficients corresponding to n and i, that

is b i ≡ n! i!(n -i)! .
The program gaussla.m, described in appendix, is a Matlab function that gives as output the values of the zeros x i and weights w i given n as the unique input. It used the standard roots.m file that calculate the roots of a polynomial using eigenvalue method. It also used the factorial function that calculates the factorial value of a given number. For example with n = 4, [x,w]=gaussla(4) returns : and where w i are some weights defined by

w = 9.39507091230114 x = 0.
w i = 2 n+1 (n!) √ π (H n (x i )) 2 .
Again, after some algebras it turns out that the polynomial H n (x) can be rewritten by :

H n (x) = (-1) n e x 2 d n dx n (e -x 2 ) (3.2) = (-1) n ñ i=1 b i x n-2(i-1)
where ñ ≡ 1 + n 2 if n is an even number, and ñ ≡ n+1 2 if n is an odd number. The b i are more complicated to find and require some specific procedures. Indeed, the coefficients b i are set as follows :

n 0 1 1 -2 2 4 - 2
(x) = -b 1 x = 2x. With n = 2, we have ñ = 2 and H 2 (x) = (-1) 2 (b 1 x 2 + b 2 ) = 4x 2 -2.
At a first glance, an explicit computational relationships between each raw and column of the above values does not seem to exist. By reversing the order of the columns at each raw, we get :

n 0 1 1 -2 2 -2 4

Intermediate column values

We use a forward recursive to compute each intermediate column values. Notice that the number of columns is given by ñ. Hence, the first new column (second column) appears at ñ = 2, that is n = 2 or n = 3.

n even

The rule to calculate the b(n, 2), with ñ ≥ 3 is as follows. First one has to compute the division value:

d 2 = |b(2, 2)/b(2, 1)| = |4/ -2| = 2
The successive values for the column 2 are given by the rule :

b(n, 2) = (-1)b(n, 1)(l n d 2 -2)
For example, with n = 4 we have l n = 3 and b(4, 2) = (-1)b(4, 1)(3 * 2 -2) = -48. Then, with n = 6 and l n = 4, we have b(6, 2) = (-1)b(6, 1)(4 * 2 -2) = 720. And so on. The same procedure has to be repeated at each new column. That is, calculate first

d n = |b(n, l n )/b(n, l n -1)|
with n = 3 and with n = 2n -2. Then, for highest value of n, use the forward recursive rule :

b(n, n ) = (-1)b(n, n -1)((l n -1)d n -d n )
For example, value of the fourth column is calculated as follows. For example, values of the third column (i.e. n = 3 and n = 5) are given by: and so on...

The program

The

Evaluation in Matlab

The evaluation of one specific function f (x) may now easily be done by creating appropriate Matlab function that calls as subroutine the gaussla.m or gausshe.m files. Hence, for example a typical matlab function, let call this file numintla.m, that evaluates an integral using Gauss-Laguerre quadrature will be defined as follows :

function [sol]=numintla(fun,nb) [xx,w]=gaussla(nb); % Call the gaussla.m file sx=size(xx,1); sol=0; for i=1:sx, x=xx(i); fx=eval(fun); sol=sol+w(i)*fx; end; sol is the value of the evaluation, n the number of points used and fun the function f (x) to be evaluated (fun must be a string). For example, to evaluate ∞ 0 e -x sin(x)dx using 6 points (n = 6), just call: sol=numintla(śin(x) ,6). The returned value is 0.50004947479768 while the exact value is 0.5.

2

  Gauss-Laguerre quadrature Evaluation of integral of the form ∞ 0 e -x f (x)dx using Gauss-Laguerre quadrature involves to calculate n i=1 w i f (x i ), where x i are zeros of the associated polynomial L n (x) = e x d n dx n (e -x x n )

e -x 2 f

 2 (x)dx using Gauss-Hermite quadrature involves to calculate n i=1 w i f (x i ), where x i are zeros of the associated polynomial H n (x) = (-1) n e x 2 d n dx n (e -x 2 )

  Last coefficient of each raw: b(n, l n ) n > 0 and evenThe last coefficient of each raw b(n, l n ), where l n is the ñ value corresponding to n, is given by b(n, l n ) ≡ 4 ln-1 . For example, with n = 4, we have l n = ñ = 3 and b(4, 3) = 4 2 = 16.n > 0 and oddWe have b(n, l n ) ≡ (-2)4 ln-1 . For example, with n = 7, we have l n = ñ = 4 and b(7, 4) = (-2)4 3 = -128.

First

  coefficient of each raw: b(n, 1) n > 2 and even Given the initial value b(2, 1) = -2, we have the following rule that can be iterated forward : b(n, 1) = (-1) n (2 + 4(l n -2))|b(n -2, 1)|. For example, with n = 8, and so l 8 = 5, knowing that b(6, 1) = -120 we have b(8, 1) = (-1) 8 (2 + 4(5 -2))|b(6, 1)| = 14 * 120 = 1680. n > 1 and odd Given the initial value b(1, 1) = -2, we have the following rule that can be iterated forward : b(n, 1) = (-1) n (2 + 4(l n -1))|b(n -2, 1)|. For example, with n = 3 we have l 3 = 2 and b(3, 1) = (-1) 3 (2 + 4(2 -1))|b(1, 1)| = -6(-2) = 12. Then with n = 5 we get b(5, 1) = (-1) 5 (2 + 4(3 -1))|b(3, 1)| = -10(12) = -120, and so on.

  First we set d 4 = |b(6, 4)/b(6, 3)| = 2/15 and then calculate b(8, 4) = (-1)b(8, 3)(3 * (2/15) -2/15) = -1 * (13440) * (8/15) = -7168. For the fifth column we calculate d 6 = 1/14 and we can obtain b(10, 6) = (-1)b(10, 5) * (3(1/14) -1/14) = -1 * (161280) * (2/14) = -23040. n odd The procedure is the same: calculate first the division value of the new column d n = |b(n, l n )/b(n, l n -1)| (with n = 2n -1), then for highest value of n use the forward recursive rule : b(n, n ) = (-1)b(n, n -1)((l n -1)d n -d n ).

d 3 =

 3 |b(5, 3)/b(5, 2)| = 0.2 b(7, 3) = (-1)b(7, 2)(3 * (0.2) -0.2)) = -1 * (-3360) * 0.4 = 1344. b(9, 3) = (-1)b(9, 2)(4 * (0.2) -0.2)) = -1 * (80640) * 0.6 = -48384.

  program gausshe.m, described in appendix, is also a Matlab function that gives as output the values of the zeros x i and weights w i given n as the unique input. With n = 4,

	[x,w]=gausshe(4) returns :
	w = 1.65068012388578 x = 0.08131283544724
	-1.65068012388578	0.08131283544724
	0.52464762327529	0.80491409000551
	-0.52464762327529	0.80491409000551


Thomas Vallée
email: thomas.vallee@univ-nantes.fr

Comparison of different Stackelberg solutions in a deterministic dynamic pollution control: the time inconsistency problem revisited

Keywords: Time inconsistency, Dynamic Stackelberg game, Pollution control JEL Codes: C7

In this paper, we study the two well-known supposed consequences of time inconsistency: the suboptimality of the time consistent solution, and the assuming increase of the follower's cost. To achieve such a goal, we study different dynamic Stackelberg solutions within a pollution control problem framework. This study is made under the assumption of different information structures, mainly we assume open-loop, feedback and closed-loop structures of information. Some of the numerical results may appear counterintuitive. Hence, there may exist some situations where a time consistent solution is optimal in comparison of the time inconsistent one. Moreover, the perfect discretionary solution may be beneficial to everyone.

Introduction

When a firm pollutes while producing, it is known that this flow of pollution will negatively affect other economic agents. If the firm is not liable to directly compensate these agents for the nuisances it causes, the production and the associated pollution levels optimal for the firm will not be optimal for society as a whole. One of the main problems in environmental economics is to find ways for a regulator to force such a firm to make socially optimal decisions, for example, through a proper use of taxes.

The problem has been extensively treated for the static case (see for example [START_REF] Pethig | Conflicts and Cooperation in Managing Environmental Resources[END_REF]). However, regulatory taxes have both short and long term consequences on the social welfare and on the firm's behavior. Taking these properly into account makes an explicitly dynamic analysis imperative. As noted by Batabyal [START_REF]Consisteny and optimality in a dynamic game of pollution control II: Monopoly[END_REF], among others, a natural way to conduct such an analysis is to model the interaction between the regulator and the firm as a dynamic Stackelberg game with the regulator as the leader.

Depending on the information structure many dynamic Stackelberg solutions do exist. In this paper, using a discrete time dynamic model of pollution control, we derive three of them, that is the open-loop, feedback and global (closed-loop) Stackelberg solutions and compare them. As simple as may seem the model, the derivation of the different dynamic Stackelberg solutions are not straightforward.

It is well-know, since the seminal works of Kydland and Prescott [START_REF] Kydland | Rules rather than discretion: The inconsistency of optimal plans[END_REF], and Barro and Gordon [START_REF] Barro | A positive theory of monetary policy in a natural rate model[END_REF][START_REF] Barro | Rules, discretion and reputation in a model of monetary policy[END_REF], that open-loop Stackelberg solutions are time inconsistent. From this literature, two others conclusions have been generally admitted. First, the discretionary solution is worst for the follower than the open-loop one with commitment. Second, the time consistent solution is suboptimal. Using numerical simulations, we show that those 1 two conclusions do not hold.

The plan of the paper is as follows. In the next section we define the pollution control model. Then in section 3, we derive the different dynamic Stackelberg solutions depending on the information's structure facing each player in the following order: first the open-loop one, second the feedback one, and third the global Stackelberg solution (that is a closed-loop structure of information). In section 4, using two numerical simulations, we compare these solutions. Finally we conclude.

The pollution control 2.1 The general model

We consider a discrete time version of the continuous time model of pollution suggested by Batabyal [START_REF]Consisteny and optimality in a dynamic game of pollution control II: Monopoly[END_REF]. There are two players: the regulator (the leader, R) and a monopolist (the follower, F ). The planning horizon is T periods, with T ≤ 20. There is no discounting. The goal of the monopolist is to maximize its cumulated profits over the T periods with respect to its choice of output. In each period t, the monopolist's revenue is given by P (q t )q t , where q t is its output in period t, and where P (q t ) is the inverse demand curves it faces.

Following Batabyal [START_REF]Consisteny and optimality in a dynamic game of pollution control II: Monopoly[END_REF], the monopolist is facing three kinds of costs associated with q t . First, a production cost wq t that is assumed to be proportional to the output. Second, the tax paid to the regulator τ t q t . And third, a cost c(x t )q t that depends on the current stock of pollution, x t . This last cost reflects the fact that the production efficiency decreases as the environment becomes more polluted. It may be or not internalized by the firm. The monopolist's optimization problem is thus given by

J F = T t=1 P (q t )q t -wq t -τ t q t -c(x t )q t → {q t } t∈[1,T ] max (2.1.1)
We assume that P ′ (q t ) < 0 and P ′′ (q t ) ≥ 0, and that c ′ (x t ) > 0, c ′′ (x t ) < 0 and c(0) = 0. Furthermore, we assume w > 0.

The regulator attempts to maximize, through its choice of tax rates, its cumulated payoff. Again, following Batabyal [START_REF]Consisteny and optimality in a dynamic game of pollution control II: Monopoly[END_REF] this payoff depends on three components. First, a function B(q t ) that represents a social benefit when tithe firm produces at the level q t . Second, a function D(x t ) which measures the damage from pollution. And finally the amount of money given by the tax τ t q t . So, the cumulated regulator's payoff is

J R = T t=1 B(q t ) + τ t q t -D(x t ) (2.1.2)
We assume that B(.) and D(.) are respectively at least C 2 and C 1 functions. Furthermore, [B ′ (q t ) > 0, B ′′ (q T ) < 0, D ′ (x t ) > 0 and D ′′ (x t ) > 0, that is the social costs of pollution are increasing in the pollution stock at an ever increasing rate. The strict concavity of B(q t ) + τ t q t is needed in order to insure the existence and uniqueness of a solution.

Finally, we suppose that x t evolves according to

x t+1 = f (q t , x t ) (2.1.3)
with x 1 given, and where f (q t ) is a differentiable function, with f ′ (q t ) > 0 and f ′′ (q t ) > 0. We also have f ′ (x t ) > 0 and f ′′ (x t ) > 0. Hence, the pollution stock in t + 1 is increasing in the pollution stock and in the firm's output in t.

For the purpose of the paper, we more specifically assume 1 :

P (q t ) ≡ a -bq t , (2.1.4) c(x t ) ≡ αx t , (2.1.5) B(q t ) ≡ γq t - q 2 t 2 , (2.1.6) D(x t ) ≡ δx 2 t 2 , (2.1.7)
x t+1 ≡ βq t + βx t .

(2. 1.8) where the coefficients a, b, α, γ, β and β are supposed to be strictly positive and with β < 1 and β < 1. The functional forms, as well as the hypotheses made earlier on the different derivatives, are standard in economic theory and will not be further justified here. The assumption β < 1 captures the fact that there is a natural resorption of the current pollution stock, at the rate (1 -β).

We may now derive the different dynamic Stackelberg solutions.

The different solutions

We assume that there is no uncertainty and that the regulator knows perfectly the different parameters of the monopolist's profits, even his cost. Furthermore, the regulator, our leader, is strong enough to force the monopolist to take as given the level of taxation.

The open-loop Stackelberg solution

This solution was first introduced by Simaan and Cruz [START_REF] Simaan | On the Stackelberg strategy in nonzero-sum games[END_REF][START_REF] Simaan | Additional aspects of the Stackelberg strategy in nonzero-sum games[END_REF] (for a more detail on it, see Başar and Olsder [START_REF] Başar | Dynamic Noncooperative Game Theory[END_REF]). To achieve the solution, the following steps are required. First, to any fixed action of the leader τ t , the reaction function of the follower is derived by maximizing the firm's payoff under the state constraint (2.1.3). Then, integrating this reaction function into the leader's payoff and minimizing again under the state constraint, gives the optimal action of the leader which induces an optimal action for the follower. As noticed by Simaan and Cruz [START_REF] Simaan | Additional aspects of the Stackelberg strategy in nonzero-sum games[END_REF], latter by Kydland [START_REF] Kydland | Noncooperative and dominant player solutions in discrete dynamic games[END_REF] and popularized by Kydland and Prescott [START_REF] Kydland | Rules rather than discretion: The inconsistency of optimal plans[END_REF], this solution is time inconsistent. (abq t )q twq tτ t q tαx t q t (3.1.1) subject to

x t+1 =βq t + βx t (3.1.2)
Let define the firm's Hamiltonian by

H F (q t , x t , p F t+1 ) ≡ J F t + p F t+1 (βq t + βx t (3.1.3)
By using the first order conditions required to maximize this Hamiltonian function, and after some algebras, we get

q = a -w -τ t -αx t + βp F t+1 2b (3.1.4) x t+1 = β(a -w -τ t -αx t + βp F t+1 ) 2b + βx t (3.1.5) p F t+1 = -α(a -w -τ t -αx t ) 2b + ( β - αβ 2b )p F t+1 (3.1.6)
1 Some others specifications are possible, see Batabyal [START_REF]Consistency and optimality in a dynamic game of pollution control I: Competition[END_REF][START_REF]Consisteny and optimality in a dynamic game of pollution control II: Monopoly[END_REF] 3 with initial and final condition p F T +1 = 0 and x 1 given. The stock of pollution at the period T + 1, x T +1 is free. One reason to let it free is that the regulator may not know what is or not an acceptable final level of pollution.

This above set of equations defines the reaction function of the monopolist (follower) to any announced tax path. Replacing (3.1.4) into J L t , and given (3.1.5) and (3.1.6), we may now solve the regulator's problem defined by the following Hamiltonian

H R (τ t , p L t+1 , p F t+1 , x t , µ t ) ≡ (γ + τ t )(a -w -τ t -αx t + βp F t+1 ) 2b - 1 2 ( a -w -τ t -αx t + βp F t+1 2b ) 2 - δx 2 t 2 + p R t+1 ( β(a -w -τ t -αx t + βp F t+1 ) 2b + βx t ) + µ t ( -α(a -w -τ t -αx t ) 2b + ( β - αβ 2b )p F t+1 ) (3.1.7)
Then we know from Başar and Olsder [START_REF] Başar | Dynamic Noncooperative Game Theory[END_REF] that the open-loop Stackelberg solution is given by the resolution of the following first-order conditions:

∂H R t ∂τ t = -γ -τ t -βp R t+1 + αµ t 2b + (a -w -τ t -αx t + βp F t+1 )(1 + 2b) 4b 2 = 0 (3.1.8) x t+1 = ∂H R t ∂p L t+1 = β(a -w -τ t -αx t + βp F t+1 ) 2b + βx t (3.1.9) p R t = ∂H R t ∂x t = p R t+1 (2b β -αβ) -α(γ + τ t ) + α 2 µ t 2b + α(a -w -τ t -αx t + βp F t+1 ) 4b 2 -δx t , (3.1.10)
p F t = ∂H R t ∂µ t = -α(a -w -τ t -αx t ) 2b + ( β - αβ 2b )p F t+1 (3.1.11) µ t+1 = ∂H R t ∂p F t+1 = (γ + τ t )β + (2b β -αβ)µ t 2b + β(βp R t+1 -2b(a -w -τ t -αx t + βp F t+1 )) 2b (3.1.12)
with x 0 given, and

µ 1 = 0 (3.1.13)
The boundary condition µ 1 = 0 is directly related to p F T +1 = 0. Furthermore, we have p R T +1 = 0. As known, the open-loop Stackelberg solution is time inconsistent, since a reoptimization latter in time, at period k for example, will give again to set µ k = 0 although initially calculated, at period 1, we have µ k = 0. Anyway, these above necessary conditions, after some algebras and following Medanic [START_REF] Medanic | Multilevel Stackelberg strategies in linear-quadratic systems[END_REF] give us to solve an augmented discrete Hamiltonian matrix (i.e. with a tracking matrix) of the form:

xt+1 pt = A B C A xt pt+1 + D E (3.1.14)
Where A, B, C are some 2 × 2 matrices, D is a 2 × 1 matrix and xt and pt are some 2 × 1 vectors defined by:

A ≡ β -βα 4b+1 -βα 4b+1 -βα 4b+1 β -βα 4b+1 , B ≡ β 2 4b+1 β 2 4b+1 β 2 4b+1 β 2 4b+1 , C ≡ α 2 4b+1 -δ α 2 4b+1 α 2 4b+1 α 2 4b+1 , D ≡ β(a-w+γ) 4b+1 β(a-w+γ) 4b+1 , E ≡ -α(a-w+γ) 4b+1 -α(a-w+γ) 4b+1 , xt ≡ x t µ t , and pt+1 ≡ p R t+1 p F t+1 .

Resolution

To solve this tracking problem defined above we use the sweep method (see Bryson and Ho [6]). That is, we assume a linear relation between the costate and the state vectors:

pk = S k xk -g k (3.1.15)
Thus, using this into the augmented Hamiltonian matrix we first get an expression for x k+1 :

xk+1 = (I 2×2 -BS k+1 ) -1 (Ax k -Bg k+1 + D) (3.1.16)
Then using (3.1.16) and (3.1.15) into the definition of p k+1 as given by the augmented Hamiltonian matrix, and equating both sides we finally get the difference equations:

S k = C + AS k+1 (I 2×2 -BS k+1 ) -1 A, (3.1
.17)

g k = AS k+1 (I 2×2 -BS k+1 ) -1 (Bg k+1 -D) + Ag k+1 -E, (3.1.18)
where the first equation is the so-called Riccati difference equation, and the second one defines a tracking difference equation.

The boundary conditions are:

The optimal discretionary open-loop Stackelberg solution

As we said, the open-loop Stackelberg solution is time inconsistent. That is, whatever is the optimal sequence of taxation {τ * } T 1 calculated at time t = 1, it will not be optimal to pursue with it at time t = 2. Rather, the regulator is induced to solve the truncated problem starting at time t = 2 and ended at t = T which give him a 5 new optimal sequence of taxation {τ * * } T 2 . Again at t = 3 this sequence will be suboptimal, and so on.

So, while {τ * } T i is the optimal open-loop sequence of taxation for the problem starting at time t = i and ended at time t = T , let {τ * i } T i be the first component of this sequence (and also unique one for the case where i = T ). Then, the optimal discretionary sequence of taxation, realized ex post, is

{τ d * t } t∈[1,T ] = ({τ * t } 1 1 , {τ * t } 1 2 , ..., {τ * t } 1 T -1 , {τ * t } 1 T ) (3.1.22)
In the economic literature, such a discretionary policy is generally assumed to be worst for the follower regardless to the initial committed strategy that is {τ * } T 1 . But, as we will see, such a conclusion may not hold. In consequence, both players, monopolist and regulator, may gain by using such a discretionary policy. If so, the monopolist may rationally accept ex post any revisions of the initial strategy.

The feedback solution

To solve the game, under the feedback structure of information assumption, we use the dynamic programming method with appropriate value functions (see Başar and Olsder for more details [START_REF] Başar | Dynamic Noncooperative Game Theory[END_REF]). Recall that this solution is time consistent by construction.

Let T be the last period of the problem. Since the level of pollution x T +1 is free, the reaction function of the monopolist at the last period is directly given by the resolution of arg max

qT ∈ℜ J F T (3.2.1) That is q * T = a -w -τ T -αx T 2b (3.2.2)
where x T is a known fixed value. Then the problem facing the regulator is simply given by: arg max

τT ∈ℜ B(q * T ) + τ T q * T -D(x T ) (3.2.3)
where q * T is given by (3.2.2). The maximum is obtained when

τ * T = (1 + 2b)(a -w -αx t ) -2bγ 1 + 4b (3.2.4)
For the problem starting at the period T -1, the value functions are defined by

V F (T -1, T ) = [arg max qT -1 J F T -1 ] + J F * T , (3.2.5) V R (T -1, T ) = [arg max τT -1 J R T -1 ] + J R * T , (3.2.6)
where J F T and J R T are known and supposed to be defined in a linear quadratic form

J L * T = P T x 2 T + p T x T + n T , (3.2
.7a)

J F * T = PT x 2 T + pT x T + n T (3.2.7b) Using the definition x T = βq T -1 + βx T -1 into J F * T
and maximize the value function for any fixed τ T -1 gives an optimal action for the monopolist for the period T -1. Integrating this into the value function of the regulator, using again the state equation definition and maximizing over all possible τ T -1 , we can write the results in some general specific forms

τ t = K t x t + k t , (3.2.8
)

q t = Kt x t + kt , (3.2.9)
x t+1 = Ω t x t + β kt , (3.2.10)

J R t = P t x 2 t + p t x t + n t , (3.2.11)
J F t = Pt x 2 t + pt x t + ñt , (3.2.12)
for all t ∈ [1, T ] and with

K t = -α + 2α(β 2 (P t+1 + Pt+1 ) -b) + 2β β( Pt+1 -2b(P t+1 -Pt+1 ) -2β 2 P 2 t+1 ) 1 + 4b -2β 2 (P t+1 + 2 Pt+1 ) , (3.2.13)
k t = -a -2ab + 2bγ + 2bβp t+1 + 2aβ 2 P t+1 -β pt+1 (1 + 2b -2β 2 P t+1 ) -1 -4b + 2β 2 (P t+1 + 2 Pt+1 ) + 2β 2 Pt+1 (a -γ -w -βp t+1 + β pt+1 ) + w -2bw -2β 2 wP t+1 -1 -4b + 2β 2 (P t+1 + 2 Pt+1 (3.2.14) Kt = α + 2β β(P t+1 + Pt+1 ) 1 + 4b -2β 2 (P t+1 + 2 Pt+1 ) , (3.2.15) kt = a + γ -w + β(p t+1 -pt+1 ) 1 + 4b -2β 2 (P t+1 + 2 Pt+1 ) , (3.2.16
)

Ω t = β Kt + β, ( 3
.2.17)

P t = -K2 t 2 + K t Kt - δ 2 , (3.2
.18)

p t = γ Kt -Kt kt + K t kt + Kt k t , (3.2.19
)

n t = -k2 t 2 + γ kt + k t kt , (3.2
.20) Pt = -b K2 t -α Kt -K t Kt , (3.2
.21) pt = a Kt -2b Kt kt -w Kt -α kt -K t kt -Kt k t , (3.2
k T = (1 + 2b)(a -w) -2bγ 1 + 4b , (3.2.25) KT = -α 1 + 4b , (3.2.26) kT = a -w + γ 1 + 4b . ( 3

.2.27)

To get the optimal feedback Stackelberg solutions, one must first solve off-line the set of equations (3.2.13)-(3.2.23) using the terminal conditions, and then compute on line the values of τ t , q t and x t .

The global Stackelberg solution

Assume that the structure of the information facing each player is a closed-loop one. That is, the leader has a perfect knowledge of all the past and current values of the state and controls which lead him to also know the monopolist's actions. In such an information structure, the regulator may try to find an incentive strategy such that he can reach his global optimum (i.e. optimum optimorum).

This optimum optimorum is assumed to be unique. Let call (q * t , τ * t ), ∀t ∈ [1, T ], the pair of actions that globally maximized J R t . Following Başar and Olsder [START_REF] Başar | Dynamic Noncooperative Game Theory[END_REF], we know that this pair of actions is directly given by the resolution of max

{τ } T 1 ,{q} T 1 J R ({q} T 1 , {τ } T 1 ) (3.3.1)
A necessary condition to solve (3.3.1), by using the first order conditions, is the strict concavity of J R t (q t , τ t ) in q t and τ t , ∀t. That is we need no existence of singularity 2 . Pity, J R t (q t , τ t ) is singular in τ t , ∀t. Hence, a direct optimization is not possible.

One way to avoid this problem, is to add a specific constraint on τ t or q t in order to reintroduce τ t in the maximization's problem. Such a constraint may be the willingness that the monopolist's profits get close to a given level. For example, the regulator may want to reduce these profits to zero 3 . To require J F t = 0 involves that either q t = 0, or

q t = a -w -αx t -τ t b . (3.3.2)
Obviously, q t = 0 must be excluded as a possible choice.

The regulator's problem may be defined by the Hamiltonian-Lagrangian function

L R t = J R t + p R t+1 (βq t + βx t ) + λ t ( a -w -αx t -τ t b -q t ) (3.3.3)
To maximize over τ t and q t this function, one must solve the following set of first order conditions

∂L R t ∂τ t = q t - λ t b = 0, (3.3.4) ∂L R t ∂q t = γ -q t + τ t + βp R t+1 -λ t = 0, (3.3.5) x t+1 = ∂L R t ∂p R t+1 = βq t + βx t , (3.3.6)
p R t = ∂L R t ∂x t = -δx t + βp R t+1 -αλ t , (3.3.7)
∂L R t ∂λ t = a -w -αx t -τ t b -q t = 0. (3.3.8)
After some algebras, we get

λ t = b(a -w + γ -αx t + βp R t+1 ) 2b + 1 . (3.3.9)
Using this, the following augmented Hamiltonian system has to be solved

x t+1 p R t = β -αβ 2b+1 β 2 2b+1 -δ + α 2 b 2b+1 β -αbβ 2b+1 x t p R t+1 + β(a-w+γ) 2b+1 -αb(a-w+γ) 2b+1 (3.3.10)
with the boundary conditions p R T +1 = 0, and x 1 given. (3.3.11) By assuming a linear relationship between the co-state and the state, p R t = K t x tg t , the following scalar 2 Namely the condition that J R T (qt, τt) is at least a C 2 function in qt and τt. 3 We chose this zero profits' constraint for simplicity purpose. Hence, the constraint J R t = 100, ∀t, involves

qt = a -τt -w -αxt --4 b n + (-a + τt + w + αxt) 2 2 b ou qt = a -τt -w -αxt + -4 b n + (-a + τt + w + αxt) 2 2 b
which lead the regulator's problem to be a non linear optimal control one hard to solve analytically.

Riccati and tracking difference equations have to be solved off line, backard in time

K t = -δ + α 2 b 2b + 1 + ( β -αbβ 2b+1 )K t+1 ( β -αβ 2b+1 ) 1 -β 2 Kt+1 2b+1 , (3.3.12) g t = αb(a -w + γ) 2b + 1 + ( β -αbβ 2b+1 )K t+1 ( β 2 gt+1-β(a-w+γ) 2b+1 ) 1 -β 2 Kt+1 2b+1 + ( β - αbβ 2b + 1 )g t+1 . (3.3.13)
with the terminal conditions

K T = -δ + α 2 b 2b + 1 , K T +1 = 0, (3.3.14) g T = αb(a -w + γ) 2b + 1 , g T +1 = 0. (3.3.15)
Once these off line values are found, the optimal sequences {x * } T 1 , {p R * } T 1 , {λ * } T 1 , {τ * } T 1 and {q * } T 1 can be calculated on line. Recall that {τ * } T 1 and {q * } T 1 achieve the optimum optimorum of the regulator given the zero-profit constraint.

Next, the problem facing the regulator is to find and announce at the beginning of the game an optimal incentive strategy such that the monopolist implements the sequence {q * } T 1 . Since the regulator know either directly q t , ∀t, or at least may calculate it from its knowledge of x t , and following Başar and Olsder [START_REF] Başar | Dynamic Noncooperative Game Theory[END_REF], we know that a candidate incentive strategy, call it θ, is

τ t ≡ θ t (q t ) = τ * t + k t (q * t -q t ) (3.3.16)
where τ * t and q * t are the desired actions from the viewpoint of the regulator, and are some known values. To proceed with θ, the regulator has to find the sequence {k} T 1 such that the monopolist cannot do better than {q * } T 1 to which the regulator responds by {τ * } T 1 . If such a sequence of incentive strategies exists, the global Stackelberg solution is time consistent by hypothesis as it reaches the optimum optimorum of the regulator. Since θ t (q t ) is a known function, the problem facing the monopolist is a standard optimal control one. Furthermore, given that there is no uncertainty, the solution will be the same whatever the information structure, open-loop or feedback, faced by the monopolist is. For simplicity purpose, we derive the solution using the dynamic programming method. So, let the incentive strategy for the last period be

θ T = τ * T + k T (q * T -q T ). (3.3.17)
At this last period, the monopolist problem involves to solve arg max qT

J F T (q T , θ T ) (3.3.18)
The first order condition is

q T = a -τ * T -k T q * T -w -αx T 2b -2k T . (3.3.19)
Recall that the equality q T = q * T is desired. Let k * T be an "incentive coefficient" such that this equality holds. Its value is given by

k * T = -(a -w -αx T -τ * T -2bq * T ) q * T . (3.3.20)
We may easily guess the sign of k * T . It should be positive since the paire (τ * T , q * T ) is calculated given a non-profit constraint which means that the monopolist, given τ * T , should not be able to produce more (i.e. q T ≥ q * T ⇒ J F T (τ * T , q T ) < 0). Since the monopolist may only decide to produce less, a lower value of q T should be associated to an increase of τ T in order to induce the monopolist to choose q * T . In consequence, k * T > 0 is required.

Anyway, the last period payoff of the monopolist can be defined as follows

J F T = (a -bq * T )q * T -wq * T -αx T q * T -τ * T q * T = PT x 2 T + pT x T + ñT , (3.3.21)
where

PT = 0, pT = -αq * T , ñT = (a -bq * T )q * T -wq * T -τ * T q * T .
One may check that θ T (k * T ) also induces the regulator to implement τ * T . The regulator payoff can be rewritten as follows

J R T = P T x 2 T + p T x T + n T (3.3.22)
with

P T = -δ 2 , p T = 0, n T = γq * T - q * 2 T 2 + τ * T q * T .
Using a similar procedure we used to get the feedback Stackelberg solution, one can get the following general forms of the closed-loop Stackelberg solution

q t = a -w -αx t + β pt+1 -τ * t -k t q * t 2b -2k t , (3.3.23)
k * t = -(a -w -αx t + β pt+1 -τ * t -2bq * t ) q * t , (3.3.24)
x t+1 = βq * t + βx t , (3.3.25)

J R t = P t x 2 t + p t x t + n t , (3.3.26)
J F t = Pt x 2 t + pt x t + ñt . (3.3.27)
where

P t = -δ 2 , p t = 0, n t = γq * t - q * 2 t 2 + τ * t q * t . Pt = 0, pt = -αq * t , ñt = (a -bq * t )q * t -wq * t -τ * t q * t .
Remark : it is possible that for some values of the parameters, we have k * t = b for some t. Then as easily seen from (3.3.19) or (3.3.23), the problem facing the monopolist becomes singular. In such a case, the optimal level of production may not be obtained by (3.3.19) or (3.3.23). In fact the optimal level of production is given by

q t =    a -w -αx t + β pt+1 -τ * t -k t q * t 2b -2k t if k * t = b, q * t if k * t = b.
(3.3.28)

4 Some numerical comparisons of the solutions

The results presented here were obtained for the following values of the parameters: a = 150, b = 5, w = 2, α = 2, δ = 3, and γ = 5.

The initial level of pollution is set to x 1 = 1. Two numerical simulations are run. In the first one, we set β = 0.4 and β = 0.5, and in tghe second one β = 0.8 and β = 0.8 are used.

4.1 First case: β = 0.4 and β = 0.5

Logically the best solution, from the regulator viewpoint, is the global one (table 1 and figure 1), and it is the worst for the monopolist since its profits reduce to zero (table1 and figure 2). This solution involves the higher levels of pollution 4 , tax and production. Recall that his global Stackelberg solution is time consistent. Quite surprising is that the time consistent feedback solution does also better than the open-loop one, with or without commitment (figure 1 and table 1). It is generally assumed that the problem of the time consistent solution is its suboptimality in respect of the discretionary one (cf. Kydland and Prescott [START_REF] Kydland | Rules rather than discretion: The inconsistency of optimal plans[END_REF], Barro and Gordon [2,[START_REF] Barro | Rules, discretion and reputation in a model of monetary policy[END_REF]). What we learn from this simple model it's that there is no way it should be always the case when the follower has a real payoff function and not a very restrictive one5 . The level of pollution is directly related to the regulator's welfare. And since all others variables are connected each others, we found the same order of the solutions in the figures. Hence, higher welfare will imply higher pollution, and so a higher price and production.

Solutions

As the global solution involves zero-profits for the monopolist, one may wonder why the monopolist will still produce something ? Obviously, the regulator may accept some profit for the monopolist by allowing a little more pollution. That is our global solution is based on a non-profit constraint. All constraints that will involve a level of pollution between this one and the one obtained under the feedback solution will still allow the global Stackelberg solution to be the first one.

The two time consistent solutions mainly differ because of the level of taxation, this level is higher with the incentive solution (global Stackelberg) since the profit must be reduce to zero.

Another important conclusion is that, in an open-loop information structure, the discretionary solution is better for everyone than to stay committed to the initial announcement (figures 1 and 2 and table 1). In such a case, we don't see any reason why this discretionary solution should involve some loss of credibility, since the monopolist may be aware that to believe in a likely recalculated sequence of taxations will get him in a better position after. Then he may optimally believe an initial sequence of taxation knowing that the regulator will not continue with it latter.

4.2 Second case: β = 0.8 and β = 0.8

The simulation provides the same kinds of comments. That is, and the more important one, the monopolist will benefit from a not-committed regulator's policy to the open-loop initial solution (table 2 and figure 4).

For the regulator, the feedback time consistent solution is no more better than the optimal discretionary one (figure 3 andtable 2). But these solutions are very closed. Finally, it seems that the gain from not staying committed to an initial open-loop solution (by using the optimal discretionary solution) is always quite small. So, the incentive to deviate is not very strong (tables 1 and 2).

Conclusion

In this paper, we derived the different possible Stackelberg solutions of a leader-follower pollution game. The different solutions are well-known, mainly because of the work of Başar and Olsder [START_REF] Başar | Dynamic Noncooperative Game Theory[END_REF]. But despite this fact, some misunderstandings still exist concerning the comparison of these solutions. We underline the incorrectness of two of them: the suboptimality of the time consistent solution, and the assuming increased cost on the follower when the leader use a discretionary policy.

Hence, with one particular numerical simulation we presented, we found that the time consistent solution is the best one for the leader. Moreover, it is possible to find a simulation such as this conclusion also holds for the follower. The gain for both players of using optimal discretionary solution was underlined. This result is closely related to the fact that a cheating-by-second play strategy may also be a good strategy for both players (see Vallée, Deissenberg and Başar [START_REF] Vallée | Optimal open loop cheating in dynamic reversed linear quadratic stackelberg games[END_REF]). Finally, we concluded on the very small advantage of using such a solution.

Of course, those results were found with a very specific dynamic game model. Another one may give opposite results. Some more theoretical understandings of the different dynamic solutions are needed if we want, for example, to know exactly when and why a time consistent solution may be suboptimal or not. Such a project is a currently research.

A Graphics

We use the following abbreviations and notations for the graphics:

• ol (-) Open-loop solution,

• old (o) Open-loop discretionary solution,

• fd (+) Feedback solution,

• cl (--) Closed-loop solution (myopic and nonmyopic cases),

  Let the time interval be[1, T ]. To any fixed τ t , t ∈ [1, T ] the firm solves arg max qt∈ℜ * T t=1

  S T +1 = 0 2×2 , and g T +1 = 0 2×1 . (3.1.20) From the boundary conditions we get S T = C, and , g T = E. (3.1.21) and so on. Once, the computation off line, backward in time, of the different values of S k and g k are made, the values of xt and pt follows. These values automatically give us the ones of x t , µ t , p R t and p F t , for all t ∈ [1, T ]. The optimal open-loop Stackelberg actions are directly given after by (3.1.8) and (3.1.4).

Figure 2 :

 2 Figure 2: Evolution of J F t

Figure 4 :

 4 Figure 3: Evolution of J R t

A. 1 Figure 8 :A. 2 Figure 12 :

 18212 Figure 5: Pollution stock

Table 1 :

 1 Cumulated Payoffs

	J R c	J F c

The reader is implicitly refereed to the corresponding graphics that are shown in appendix.

These literatures are based on some specific Stackelberg games where the follower has a kind of cheating aversion cost function.


Roar Adland
email: roar.adland@nhh.no

Pierre Cariou
email: pierre.cariou@kedgebs.com

LEMNA Francois-Charles Wolff

Comparing transaction-based and expert-generated price indices in the market for offshore support vessels

Keywords: Hedonic pricing, fixtures, offshore, dayrates, shipbrokers, heterogeneous contracts

We propose a methodology to generate transaction-based price indices in a heterogeneous physical market and to compare such indices with expert-generated price indices. We apply our framework to the offshore services market using a unique dataset of around 8,000 transactions. We find that the deviation between our transaction-based indices and shipbrokers' market indices is positively impacted by the level of dayrates and negatively impacted by the volume of transactions. Our approach is useful for the quality assurance of expert-generated price indices in heterogeneous physical markets such as transportation services.

Introduction

In most physical markets, certainly those where transactions occur in decentralised Over-the-Counter (OTC) markets, the compilation, standardization and dissemination of market price indices is usually provided by price-reporting agencies or intermediaries. In the markets for transportation services, well-known examples are the Baltic Exchange for drybulk and tanker freight rates and price indices compiled by shipbrokers, with perhaps the best known example being the variety of freight rate indices published by Clarkson Research (2017).

In virtually all these cases, the benchmark price (that is, the freight rate or day rate) is set at regular intervals, typically at a daily/weekly frequency, by asking a shipbroker or a panel of such intermediaries to submit their subjective assessment of the prevailing price. The expert's opinion is based on both public information (such as the price agreed in recent transactions or fixtures1 ) and private information (such as knowledge of ongoing negotiations between buyers and sellers). The degree of subjectivity in generating indices is likely to rise when markets are characterised by a limited number of transactions and by some heterogeneity in fixtures characteristics (cargo size, load date, origin/destination port, vessels characteristics) as discussed in [START_REF] Alizadeh | Vessel and voyage determinants of tanker freight rates and contract times[END_REF] or [START_REF] Adland | The influence of charterers and owners on bulk shipping freight rates[END_REF][START_REF] Adland | What makes a freight market index? An empirical analysis of vessel fixtures in the offshore market[END_REF]. [START_REF] Veenstra | Price indices for ocean charter contracts[END_REF] discuss the issues inherent in such expert-generated price indices as "there is very little consensus on what type of information is required by practitioners" and "the relationship between the available data set on market transactions (fixtures) and the economic indicators is a black box". It is also well known that publicly available information only covers a limited part of overall market activity (Veenstraa and van [START_REF] Veenstra | Price indices for ocean charter contracts[END_REF] and that some of these transactions may relate to the speculative re-letting of ships and not end-user demand [START_REF] Laulajainen | The Geographical Foundation of Dry Bulk Shipping[END_REF]. Despite these potential drawbacks, empirical research relies heavily on brokers' freight market indices throughout the literature, for instance in continuous-time models [START_REF] Bjerksund | Contingent claims evaluation for mean-reverting cash flows in shipping[END_REF][START_REF] Tvedt | Valuation of VLCCs under income uncertainty[END_REF][START_REF] Adland | The non-linear dynamics of spot freight rates in tanker markets[END_REF][START_REF] Koekebakker | Are spot freight rates stationary?[END_REF][START_REF] Adland | A discrete-time stochastic partial equilibrium model of the spot freight market[END_REF][START_REF] Adland | Price dynamics in the market for Liquid Petroleum Gas transport[END_REF][START_REF] Poblacion | The stochastic seasonal behavior of freight rate dynamics[END_REF]Poblacion, , 2017) ) or timeseries models (Kavussanos, 1996;[START_REF] Berg-Andreassen | Some properties of international maritime statistics[END_REF]Frances and Veenstra, 1997;Kavussanos and Alizadeh, 2001). Our study also complements the very limited literature on the market for Offshore Supply Vessel (OSV) services [START_REF] Aas | The role of supply vessels in offshore logistics[END_REF]Kaiser, 2014;Eskandari and Mahamoodi, 2016) and, in particular, former work by [START_REF] Adland | What makes a freight market index? An empirical analysis of vessel fixtures in the offshore market[END_REF] on the identification of determinants of OSV price indices.

Reliable price indices of heterogeneous goods, such as real estate (Jiang et al., 2014) or fine art [START_REF] Bocart | Volatility of price indices for heterogeneous goods with applications to the fine art market[END_REF], are important for two main reasons. First, brokers reporting their price assessments may collude to manipulate the price index higher or lower than the true market price. This has led the International Organisation of Securities Commissions (IOSCO) to introduce principles to emphasise a focus on actual transactional data, although it accepts that in certain circumstances an expert view may be a valid method to provide a benchmark (Baltic Exchange, 2014). Second, in the freight markets, the increasing use of digital auctioning platforms for the chartering of vessels means that buyers and sellers are by-passing brokers and entering into transactions directly2 . If such platforms become successful, then a larger share of transactions will be outside the control of traditional brokers, reducing their ability to generate representative indices.

The contribution of our paper is twofold. From a methodological viewpoint, we use a hedonic pricing framework to show how a shipping index based on transactional data can be built so that it provides an informational content equivalent to that produced by a shipbroker index. From an empirical perspective, we investigate how our transaction-based indices behave compared to shipbrokers' indices using about 8,000 fixtures available for the North Sea spot offshore market from 2000 to 2016. We focus on this market because it possesses an interesting combination of high volatility, high liquidity and fleet heterogeneity, which may impact the use of judgments by experts. This offers the possibility to assess the validity of our approach and to test potential determinants (like market liquidity and volatility) that may explain the gap between pure transactional-indices and broker-indices.

The remainder of the paper is organized as follows. In Section 2, we present the hedonic pricing framework to generate different market price indices based on transaction data. Section 3 discusses the data available to construct a representative transactional market index for the OSV North Sea market and presents the resulting price indices derived from our hedonic pricing framework. In Section 4, we discuss how our indices perform both in comparison to each other and to indices generated by experts (shipbrokers), and suggest potential reasons for the observed differences. Section 5 concludes.

Estimation of broker-equivalent market indices

Two main alternative methods exist to generate price indices from heterogeneous transactions: the hedonic regression method and the repeat sales method (Jiang et al., 2014).

The hedonic approach [START_REF] Rosen | Hedonic prices and implicit markets: Product differentiation in pure competition[END_REF]) -which we follow here -is to regress the price of each good on its own characteristics in order to control for variations due to observable differences between heterogeneous goods [START_REF] Bocart | Volatility of price indices for heterogeneous goods with applications to the fine art market[END_REF].

We propose a methodology to construct a market index expected to be similar in its informational content to that provided by a shipbroker index. Let us consider a period of 𝑇𝑇 months with 𝑡𝑡 = 1, … , 𝑇𝑇 and let 𝕊𝕊 𝑡𝑡 be the value of the broker index at month 𝑡𝑡. The available information set from the broker is 𝕊𝕊 = {𝕊𝕊 1 , … , 𝕊𝕊 𝑇𝑇 }2 F 3 , where each value 𝕊𝕊 𝑡𝑡 represents the broker's estimate of the market freight rate that is published at date 𝑡𝑡. For transaction data corresponding to fixtures between ship owners and charterers, we denote each fixture by ℱ 𝑖𝑖 for the 𝑁𝑁 transactions over the period of length 𝑇𝑇. The 𝑖𝑖 th fixture is signed at time 𝑡𝑡 so that ℱ 𝑖𝑖 refers in fact to ℱ 𝑖𝑖(𝑡𝑡) . For each fixture ℱ 𝑖𝑖 , we have information on the freight rate 𝑓𝑓 𝑖𝑖 , on the characteristics of the contract 𝑋𝑋 𝑖𝑖 (like vessel age, vessel capacity or contract duration) and on the charterer 𝑐𝑐3 F 4 . In what follows, we turn to a hedonic price regression à la [START_REF] Rosen | Hedonic prices and implicit markets: Product differentiation in pure competition[END_REF] The average freight rates can be generated from the following linear regression model:

𝑓𝑓 𝑖𝑖 = ∑ 𝔽𝔽 𝑡𝑡 * 𝕝𝕝 𝑡𝑡 𝑇𝑇 𝑡𝑡=1 + 𝜀𝜀 𝑖𝑖 (1)
where 𝕝𝕝 𝑡𝑡 are dummy variables such that 𝕝𝕝 𝑡𝑡 = 1 at date 𝑡𝑡 and 𝕝𝕝 𝑡𝑡 = 0 otherwise, 𝔽𝔽 𝑡𝑡 are parameters to estimate, and 𝜀𝜀 𝑖𝑖 is a residual perturbation with 𝐸𝐸(𝜀𝜀 𝑖𝑖 ) = 0 and 𝑉𝑉(𝜀𝜀 𝑖𝑖 ) = 𝜎𝜎 2 .

When the linear model (1) does not include a constant, each parameter 𝔽𝔽 𝑡𝑡 is exactly equal to the average 𝑓𝑓 ̅ 𝑡𝑡 . Our first market index is hence 𝔽𝔽 𝑎𝑎𝑎𝑎𝑎𝑎 = {𝔽𝔽 1 , … , 𝔽𝔽 𝑇𝑇 } where 𝑎𝑎𝑎𝑎𝑎𝑎 stands for average. Such market index is straightforward to implement and gives the possibility to provide some standard errors 𝜎𝜎(𝔽𝔽 𝑡𝑡 ) for each parameter 𝔽𝔽 𝑡𝑡 and therefore to construct a confidence interval. However, 𝔽𝔽 𝑎𝑎𝑎𝑎𝑎𝑎 does not consider the fact that fixtures are heterogeneous [START_REF] Adland | The influence of charterers and owners on bulk shipping freight rates[END_REF].

For instance, consider the case where charterers are selecting larger vessels compared to those hired during the previous months. As there exists a positive correlation between carrying capacity and freight rate, the average price for that month will rise. Nevertheless, this increase in the average market index is only due to a change in the vessel characteristics for this specific month. To account for this composition effect, we extend the regression model

(1) in the following way. First, we control for observable heterogeneity by adding both the vessel and contract characteristics 𝑋𝑋 𝑖𝑖 as covariates. Second, we control for unobservable heterogeneity by adding a charterer fixed effect 𝜃𝜃 𝑐𝑐 , where 𝑐𝑐 is the charterer. The fixed effect 𝜃𝜃 𝑐𝑐 picks up all the observed and unobserved characteristics of the charterers that remain constant over time. The model to estimate is:

𝑓𝑓 𝑖𝑖 = ∑ 𝔽𝔽 𝑡𝑡 * 𝕝𝕝 𝑡𝑡 𝑇𝑇 𝑡𝑡=1 + 𝑋𝑋 𝑖𝑖𝑡𝑡 𝛽𝛽 + 𝜃𝜃 𝑐𝑐 + 𝜀𝜀 𝑖𝑖 (2)
with 𝛽𝛽 a vector of coefficients to estimate. Using (2), we obtain a new index 𝔽𝔽 𝑛𝑛ℎ = {𝔽𝔽 1 , … , 𝔽𝔽 𝑇𝑇 } where 𝑛𝑛ℎ stands for net of heterogeneity. The various 𝔽𝔽 𝑡𝑡 provide time-average values of the freight rates that account for heterogeneity in the composition of both vessels and charterers. This is an important feature, especially as our aim is to generate long-term market indices for which composition effects are likely to be substantial. For instance, due to technological development, vessels are larger and more energy-efficient by the end of the sample period.

Keeping in mind that brokers are publishing their market indices on a regular basis (usually weekly or monthly) and cannot change the values retrospectively, a potential shortcoming of the 𝔽𝔽 𝑛𝑛ℎ index is that it is the result of an ex-post assessment. Specifically, we use the whole set of fixtures 𝑓𝑓 𝑖𝑖(𝑡𝑡=1,…,𝑇𝑇) to estimate the coefficients 𝔽𝔽 𝑡𝑡 . In doing so, we assume that shipbrokers publish a market index having information for all 𝑇𝑇 periods. This is certainly misleading as, at date 𝑡𝑡, the information is only available for periods from 1 to 𝑡𝑡, but not from 𝑡𝑡 + 1 to 𝑇𝑇4 F 5 . To overcome this shortcoming, we propose a third market index where each coefficient 𝔽𝔽 𝑡𝑡 is estimated from repeated regressions, using observations from 1 to 𝑡𝑡. This means that we estimate 𝑇𝑇 -1 regressions equivalent to (2), but based on different samples, such that:

𝑓𝑓 𝑖𝑖,𝑡𝑡 ̃= ∑ 𝔽𝔽 𝑡𝑡 * 𝕝𝕝 𝑡𝑡 𝑡𝑡 t𝑡=1 + 𝑋𝑋 𝑖𝑖𝑡𝑡 𝛽𝛽 + 𝜃𝜃 𝑐𝑐 + 𝜀𝜀 𝑖𝑖 (3)
where 𝑓𝑓 𝑖𝑖,𝑡𝑡 ̃ refers to fixtures over the period from 1 to 𝑡𝑡 ̃ such that 𝔽𝔽 2 is obtained from a linear regression explaining fixtures over the first two periods, 𝔽𝔽 3 is obtained from fixtures over the first three periods, and so on. The new market index is referred as In what follows, we implement the three market indices using a database of offshore supply vessel (OSV) contracts. We compare the time profiles obtained from 𝔽𝔽 𝑎𝑎𝑎𝑎𝑎𝑎 , 𝔽𝔽 𝑛𝑛ℎ and 𝔽𝔽 𝑟𝑟𝑛𝑛ℎ to assess the consequences of the underlying methodology on the market indices. Then, we compare the indices provided by shipbrokers with market indices calculated from transaction data.

𝔽𝔽 𝑟𝑟𝑛𝑛ℎ = {𝔽𝔽 1 , … ,

Construction of market indices for Offshore Supply Vessels (OSV)

The OSV dataset

There are two main types of OSVs. Platform Supply Vessels (PSV) are designed to transport supplies and equipment from onshore bases to offshore oil and gas installations (OECD, 2014). Anchor-handling Tug Supply (AHTS) vessels are designed to tow mobile offshore drilling or production units, for handling the rigs' anchors, and sometimes as supply vessels when towing activities are low. These vessels are required during the entire life cycle of an offshore oil and gas field (from exploration, drilling, offshore development, life of field to decommissioning) and represents one of the main cost components of the upstream oil and gas industry [START_REF] Aas | The role of supply vessels in offshore logistics[END_REF].

The global OSV market consists of several regional markets -notably the North Sea, West Africa, Brazil and Southeast Asia. These markets correspond to the major offshore oil producing regions and have their own characteristics due to technical, legal, fiscal and geographic differences (ICS, 2011). For instance, vessels operating in the North Sea must be built to higher technical specifications to cope with the harsh weather conditions, with higher average dayrates than the simpler vessels required in Southeast Asia. The type, complexity and size of the rigs serviced, as well as local taxation, manning and safety regimes, will also affect the dayrates prevalent in the markets. For these reasons, there is only limited fungibility of vessels across geographical markets. The contract duration can also differ substantially, with a portion of vessels being hired for short duration (less than a month) on the spot market and the remainder for longer durations in the term market. The differences in the regional market structure can be illustrated by the descriptive statistics reported in Table 1, based on contract data provided by ODS Petrodata from 2000 to 20156 .

Insert Table 1 Out of the 25,844 transactions reported over the period, 75 per cent of contracts are on the spot market, varying from 78 per cent for AHTS to 72 per cent for PSVs. With 20,682 transactions (or 80 per cent of all transactions reported worldwide), the North Sea market is by far the most active market, followed by Southeast Asia (5.7 per cent) and USA-Gulf of Mexico (4.4 per cent). While the term market (the chartering market for contracts longer than 30 days duration) is larger than the spot market when measured by vessel employment (total ship-days per annum), the number of transactions per region is substantially lower and heterogeneity in terms of contract duration is larger, making it less suitable for our empirical approach. Similarly, expert brokers tasked with generating dayrate indices are also dependent on having information on a sufficient number of transactions to produce reliable estimates.

Accordingly, the dominance of the North Sea spot market within the global OSV market (by number of transactions) and low liquidity elsewhere means that it is the most appropriate target for our empirical analysis and comparison between transaction-based and expertgenerated indices7 .

Each fixture includes information related to the transaction, such as the reported date, the dayrate (in GBP/day), contract duration, operating region, vessel owner and charterer name, as well as the vessel's scope of work (activity). Furthermore, in the ODS Petrodata database, this commercial information has been matched with technical specifications of the vessel regarding carrying capacity (PSV deck area and deadweight), size of the engine installation (AHTS brake horsepower and bollard pull), type of propulsion (diesel electric or conventional mechanical), Year of build (Age), Design speed, Country of build and the sophistication of its Dynamic Positioning system. Additionally, it is noted whether the vessel has special design features such Ice class, Remotely operated vehicle support or Helicopter deck. While not ordinarily required, such features may be required for certain contracts where the vessel is supporting sophisticated marine operations in Arctic areas, for instance. Age is a running value calculated as the difference between the contract date and the date of vessel delivery.

There is a high degree of heterogeneity in transactions, both in terms of the main technical vessel specifications and the main activity provided under a contract [START_REF] Adland | What makes a freight market index? An empirical analysis of vessel fixtures in the offshore market[END_REF]. However, this heterogeneity is reduced when we focus on a specific sub-segment, such as a particular range of deck area for a PSV or engine size range for AHTS. As the purpose of our research is to compare transaction-based indices with expert-generated market indices produced by shipbrokers, we further need to narrow down our analysis by vessel size. This should ensure that our transaction-based price indices are directly comparable with those produced by expert shipbrokers. Within the North Sea spot market, Clarkson Research (2017) produces four indices: for AHTS vessel with 16,000-20,000 Brake Horse Power (BHP) and more than 20,000 BHP, respectively, and for PSVs with a deck area between 500-900 m2 and above 900 m2. From Table 1, AHTS with 16,000-20,000 BHP and more than 20,000 BHP represent 28 per cent and 32 per cent of transactions in the North Sea market, while for PSV between 500-900 m2 deck area and more than 900 m2, the respective shares of the North Sea market are 65per cent and 31 per cent. In our analysis we restrict our attention to the smaller vessel sizes, as they have a sufficient number of transactions throughout the sample period.

Figure 1 presents the AHTS and PSV average monthly freight rates and the average number of transactions for these two market segments. AHTS dayrates are higher than PSV dayrates and there are some large fluctuations in rates over time. The large increase in dayrates from 2005 to 2009 can be explained by the boom in the demand for offshore projects in the North Sea and to a shortage in the number of vessels available (OECD, 2014).

Insert Figure 1 This volatility, which is quite common in most shipping markets, can be attributed to the difficulty in adjusting supply to demand: first, in the short run, as it takes time to move a vessel from on area to another (to the extent that it is feasible at all) and second, in the long run, as it also takes time to build new vessels (Kalouptsidi, 2014). Another important factor explaining this volatility is the sensitivity of rig activity to oil prices [START_REF] Ringlund | Does oilrig activity react to oil price changes? An empirical investigation[END_REF] as well as the sensitivity in demand to weather conditions, which is quite unpredictable in the North Sea area [START_REF] Aas | The role of supply vessels in offshore logistics[END_REF]. The volatility in the average price can also be explained by the limited number of transactions reported during a given month. For instance, over the entire period and per month, there are never more than 40 AHTS transactions or more than 60

PSV transactions.

The reduced heterogeneity in vessel specifications when we focus on a smaller subsegment means that many of the covariates show insufficient variation such that they can be dropped from the analysis. To select the most relevant variables amongst the total 14 variables available for the AHTS or PSV market, we rely on market reports (ICS, 2011; OECD, 2014), information from brokers and operators, and the statistical properties of the variables. For instance, we use Brake Horsepower in lieu of Bollard pull for AHTS and Deck area in lieu of Deadweight carrying capacity for PSVs as these measures of size are highly correlated. We drop design speed both because of limited variability and because it need not be a good proxy for the speed at which the vessels are actually operated. We retain ice class and propulsion type as OSV operators maintain that these are important.

In our final set of covariates we are left with six variables for the AHTS market (brake horsepower; age; Remotely Operated Vehicle (ROV) support; ice class; conventional diesel; activity) and six variables for the PSV market (deck area in squared meters; age; dynamic positioning 2; ice class; conventional diesel; activity). For AHTS vessels, brake horsepower, ROV support, ice class and conventional diesel have an expected positive effect on day rates, while for PSVs, deck area, dynamic positioning, ice class, and conventional diesel are also likely to lead to a premium. Age is expected to impact dayrates negatively as the vessel gradually becomes less competitive compared to younger peers.

Table 2 presents some descriptive statistics from the ODS Petrodata dataset for the key variables. The average AHTS Brake horsepower is 17,584 BPH. The average PSV deck area is 709 m². The AHTS fleet is younger than the PSV fleet on average (5.1 years versus 10.3).

Less than 4 per cent of AHTS vessels have ROV support and 68 per cent are ice class vessels.

For PSVs, 55 per cent are classified as having a sophisticated dynamic positioning (DP2) system, 9.8 per cent have ROV support and less than 6 per cent are ice-classed vessels.

Finally, PSVs are primarily hired to perform cargo runs (73 per cent) and other supply duties (25 per cent), while AHTS are primarily used to move rigs (66 per cent) and secondarily cargo runs (15 per cent). The average contract duration is a mere 3.3 days for an AHTS and 3.1 days for a PSV, illustrating the extremely short-term nature of this spot market.

Insert Table 2 Other sources of heterogeneity can be attributed to differences in chartering practices and preferences among the largest buyers and to the limited number of vessel available on a given date, though we note that because of the very short contract duration there will be a fast turnover of vessels in the market. Nevertheless, the number of unique vessels chartered during a month is lower than the number of transactions, which further reduces the effective heterogeneity in our sample. On average, there are less than 15 distinct AHTS vessels chartered in a month versus around 40 transactions and less than 30 distinct PSVs chartered versus 60 transactions. The ten largest charterers in the North Sea AHTS 16,000-20,000 BHP segment represent 1,428 transactions (46 per cent of all transactions). For the small PSV (500-900 m2 deck area) segment, the ten largest charterers account for 2,537 transactions (48 per cent of all fixtures from 2000 to 2015) 8 . For OSV operators, there is less concentration and, over the sample period, various mergers and acquisitions mean that it is more difficult to track who operates the vessel over time9 . This led us to include a fixed effect related to charterers only.

The observations on the existence of a degree of heterogeneity, even when focusing on specific market segments, highlight the possible differences between our objective transaction-based index generation and how a human expert will subjectively evaluate market conditions. A broker may ignore certain transactions because they are not considered to be representative of the market situation or because they reflect a special relationship between the charterer and owner (filtering data). The broker may also have knowledge about offmarket transactions, ongoing negotiations or variables not included in our ODS Petrodata dataset such as the port of origin10 . While such private information or filters are not observable to researchers, these factors mean that the information set used by shipbrokers is different from the fixture dataset on which our estimations are based.

Estimation of the transactional-market based indices

In Table 3, we report estimates from the hedonic price regressions obtained for AHTS and PSV, respectively. In Model 1, we only account for the time fixed effects and obtain the set of average freight rates at the monthly level, 𝔽𝔽 𝑎𝑎𝑎𝑎𝑎𝑎 (equation 1). In both markets, the high R 2 (0.73 for AHTS and 0.82 for PSV, respectively) confirms that a large share of the variation in freight rates is explained by time fixed effects, that is, the supply and demand conditions at the time of the transaction.

Insert Table 3 In Model 2, we add the vessel and contract characteristics as well as charterer fixed effects to obtain 𝔽𝔽 nh (equation 2). The inclusion of the additional covariates leads to an increase in the R 2 by 7.1 percentage points for the AHTS market (from 0.730 to 0.782) and by 4.1 percentage points for the PSV market (from 0.822 to 0.856). For AHTS, a 1 per cent increase in a vessel's brake horsepower leads to a 1.02 per cent increase in the AHTS spot freight rate. For PSVs, a 1 per cent increase in deck area leads to a 0.23 per cent increase in freight rate. AHTS equipped with conventional diesel obtain a freight rate premium and dayrates are lower when AHTS vessels are hired for cargo runs, an activity which is outside their core purpose, and higher for rig moves. In the case of PSVs, dayrates are lower when vessels are used for cargo runs and higher when the vessel has ice class.

Next, we extend the framework by estimating the three price indices 𝔽𝔽 𝑎𝑎𝑎𝑎𝑎𝑎 , 𝔽𝔽 𝑛𝑛ℎ and 𝔽𝔽 𝑟𝑟𝑛𝑛ℎ proposed in Section 2. We present our results in Figure 2 for AHTS and in Figure 3 for PSV 11 . For the presentation, we normalize the market index to 100 in January 2000. For each market, we report both the price indices in level and then in difference with respect to 𝔽𝔽 𝑎𝑎𝑎𝑎𝑎𝑎 in order to shed light on the influence of the underlying econometric model on the time profile.

Insert Figure 2 Insert Figure 3 A large difference between the price indices produced by the different model specifications will indicate that accounting for the heterogeneity in technical vessel specifications and charterers is important, at least at particular points in time. For the AHTS market, we observe that in 2006 and 2007 the monthly indices associated to 𝔽𝔽 𝑛𝑛ℎ and 𝔽𝔽 𝑟𝑟𝑛𝑛ℎ are generally substantially lower than the monthly averages of dayrates. For a few months, the 11 The complete set of estimates (𝔽𝔽 𝑎𝑎𝑎𝑎𝑎𝑎 , 𝔽𝔽 𝑛𝑛ℎ and 𝔽𝔽 𝑟𝑟𝑛𝑛ℎ ) for AHTS and PSV is available upon request.

gap is as high as 800 points. Accordingly, for this market, we can conclude that controlling for vessel and charterer heterogeneity leads to substantial differences in the value of the market index. The situation is somewhat different for the PSV market. Our estimates show a tendency for 𝔽𝔽 𝑛𝑛ℎ and 𝔽𝔽 𝑟𝑟𝑛𝑛ℎ to be higher than average values with a difference of around 50 to 100 points. However, since 2010, the various indices are rather close in absolute values.

Another finding from Figures 2 and 3 is that there is little difference between our two indices that account for vessel and contract heterogeneity. This result is related to the way we account for transactions over time. While 𝔽𝔽 𝑛𝑛ℎ uses a unique regression to estimate the whole set of time effects {𝔽𝔽 1 , … , 𝔽𝔽 𝑇𝑇 }, the third model estimate the time effect 𝔽𝔽 𝑇𝑇 sequentially. In the long run, the coefficients 𝔽𝔽 𝑇𝑇 𝑛𝑛ℎ and 𝔽𝔽 𝑡𝑡 𝑟𝑟𝑛𝑛ℎ are expected to converge as 𝑡𝑡 increases since our specification estimates a unique set of coefficients for the vessel and contract characteristics over the period. As more and more information is available with time, this additional information has little influence on our results.

A comparison of transaction-based and broker-generated market indices

Next, we compare in Figure 4 and 5 our market indices based on transaction data to indices provided by shipbrokers for vessels with 16,000-20,000 BHP for AHTS and vessels with 500-899 m² deck for PSV12 . Compared to model 1 that was based on the average transaction prices for all AHTS or PSV (Figure 2 and 3), some differences exist between the average monthly dayrates and the dayrate indices reported by the broker.

Insert Figure 4 Insert Figure 5 The Clarkson monthly indices are less volatile than the simple index comprised of monthly averages of contracted dayrates. For instance, the maximum broker-reported average monthly rate is around GBP 100,000/day for an AHTS in the beginning of 2007, against a maximum at around GBP 195,000/day at the end of 2008 for the transaction-based index. The same applies to the PSV market, with a maximum of GBP 32,000/day at the end of 2006 for the broker price index, against GBP 73,000/day at the end of 2007 when using the average of transaction prices. An explanation of the lower volatility of our indices compared to shipbroker indices could be related to the discussion in Section 3.1, suggesting that some transactions might not be considered by brokers as representative, while our purely transactional-based approach considers all transactions without any filtering.

To ease the comparison, and as 𝔽𝔽 𝑛𝑛ℎ and 𝔽𝔽 𝑟𝑟𝑛𝑛ℎ are relatively close, we only report results for 𝔽𝔽 𝑟𝑟𝑛𝑛ℎ (model 3) and we calculate the difference 𝔻𝔻 𝑡𝑡 = 𝔽𝔽 𝑡𝑡 𝑟𝑟𝑛𝑛ℎ -ℂ𝑙𝑙𝑎𝑎𝑟𝑟𝑙𝑙𝑙𝑙𝑙𝑙𝑛𝑛𝑙𝑙 𝑡𝑡 and the associated confidence intervals. When the value 𝔻𝔻 𝑡𝑡 = 0 does not belong to the confidence interval for a period 𝑡𝑡, this means that our transaction index and the Clarkson index are significantly different. For the AHTS (Figure 4) segment, 𝔽𝔽 𝑟𝑟𝑛𝑛ℎ is relatively close to the index provided by brokers. The situation is similar for the PSV segment (Figure 5). This suggests that our index based on transaction data could provide a good proxy of the broker index. At the same time, some significant gaps are observed during specific periods of time. These differences usually correspond to periods with high volatility, particularly from 2006 to 2009 (Figure 1).

Owing to these observed differences, we finally investigate the potential reasons of such a gap. We argue that two main reasons may exist. First, the difference measured for a given month may be explained by the market level and by the degree of volatility of dayrates.

Our assumption is that at a time of high/low dayrates and high/low volatility, the broker might consider some transactions as not representative or may include information that is outside the transaction dataset (for instance, fixtures that are rumored but not publicly confirmed). This will lead to a gap as our indices are solely relying on reported transaction data, without filtering or considering the existence of potential outliers. Second, the difference may be explained by a limited number of transactions during a given month or occur at times when the utilization rate of vessels is particularly high/low. This may also lead brokers to disregard some datapoints, while we are not. In table 4, we report the results from a linear regression models testing for these two main hypotheses applied to the difference observed for the 192 months under consideration (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016).

Insert Table 4 The results show many similarities across the two vessel types. First, all coefficients for dayrates are positive and significant at the 1 per cent confidence level both for PSV and AHTS, while the volatility of dayrates does not significantly influence the gap between the Clarksons and our transaction-based market indices. Second, the number of transactions and the utilization rate of vessel reported during a given month are negatively related to the gap.

Therefore, during a month with many transactions, our transaction-based indices are closer to the indices reported by brokers, with the added advantage of objectively taking into account the heterogeneity in transactions. All told, our empirical evidence suggests that the deviations between transaction-based and expert-generated market price indices are greater in times of illiquid markets.

Concluding comments

In this paper, we have developed a methodology that enables us to extract market price indices from fixture data, subject to different assumptions about the information set. We provide, for the first time, a comparison of such transaction-based indices with those reported by a leading shipbroker company. A key issue is that there is no transparency regarding the type of information brokers are relying upon. In our contribution, we show that there can be substantial differences between expert-generated indices and those based on reported transactions data.

There are several competing hypotheses for the observed differences between the expert-generated indices and those derived from public fixture data. The most common one is that the expert shipbrokers have private information about "off market" transactions that are not in the public domain (and therefore not in our dataset). In traditional spot freight markets, such as those for crude oil tankers, this would be a valid criticism, as the share of public fixtures compared to overall fleet activity is known to be very low. However, the North Sea

OSV markets that we deal with here are highly transparent in comparison. Another possible explanation is that the differences between the indices are driven by heterogeneity in the composition of the fleet or the market players involved, or even in the composition of the experts setting the price benchmarks over time, which the expert-generated indices may not properly account for. Importantly, the proposed hedonic pricing methodology allows us to build OSV price indices that account for heterogeneity of vessel characteristics. This is crucial when dealing with an industry, such as the OSV market, where technical specifications and market composition can change substantially in the medium and long run.

We further investigated under which conditions the two indices differ. From our analysis, we identify two main situations to explain an increase in the observed gap: when dayrates are high and when the number of transactions is relatively low. We suggest that during such extreme market conditions or periods with fewer contracts, subjectivity in experts' judgement of which information to include may play a greater role.

We do not suggest that our market indices should replace those provided by shipbrokers, but that there is a need for more empirical analysis as these indicators are extremely important for both practitioners and academia. Despite the difficulties in establishing transaction-based indices in markets with limited and heterogeneous transactions, we suggest that the need for a better understanding of how price indices are generated will increase in the future for two main reasons. First, the regulatory environment requires everincreasing transparency and objectiveness of market price reporting, particularly in markets where there exist financial derivatives settled on such price indices. Second, digitalization and the possible automation of the chartering process will increase our reliance on data-driven market index generation. Most importantly, nearly all maritime economic research that deal with macroeconomic analysis (the analysis of time-series properties of freight rates) take these broker indices at face value, so that knowledge of any potential shortcomings is of crucial importance.

Finally, we note that both our transaction-based market indices and the expertgenerated price indices estimated by brokers only give an indication of the prevailing dayrate assuming that a contract can be obtained. In poor spot markets, an equally important question is whether there is actually employment for the vessels available, or whether vessels will end up unemployed for a prolonged period of time waiting for the next contract. Future work in this area should attempt to account for the demand/supply balance and develop indices that may give a better picture of the true cashflow from OSV operation after unemployment risk is accounted for.

Figure 1 .

 1 Figure 1. Average monthly freight rates and number of transactions for AHTS and PSV

Figure 2 .

 2 Figure 2. Estimated market indices AHTS

Figure 3 .

 3 Figure 3. Estimated market indices PSV

Figure 4 .

 4 Figure 4. Comparison of Clarksons and estimated market indices for AHTS

  /𝑛𝑛 𝑡𝑡 with 𝑛𝑛 𝑡𝑡 the number of fixtures recorded during time unit 𝑡𝑡.

	to
	generate several market indices.
	Our first index corresponds to the average freight rate 𝑓𝑓 ̅ 𝑡𝑡 for all fixtures occurring at
	date 𝑡𝑡 such that 𝑓𝑓 ̅ 𝑡𝑡 = ∑ 𝑓𝑓 𝑖𝑖(𝑡𝑡)

Table 1 . Description of OSV fixtures (2000-2015)

 1

	Type of OSV		Worldwide			North Sea
		Spot	Term	Total	Spot	Term	Total
	AHTS
	All	75.8%	24.2%	15,088	96.9%	3.1%	11,501
	16,000-20,000 BHP	89.11%	10.9%	3,538	97.6%	2.4%	3,186
	20,000+ BHP	95.6%	4.4%	3,873	98.1%	1.9%	3,766
	PSV
	All	71.2%	28.8%	11,736	86.3%	13.7%	9,449
	500-900 m²	76.6%	23.4%	7,043	88.7%	11.3%	5,978
	900+ m²	72.7%	27.3%	3,491	81.4%	18.6%	3,093
	Total	73.8%	26.2%	26,824	92.2%	7.8%	20,950
	Source: authors' calculations, OSV database

Table 2 . Northwest Europe OSV fixtures (2000-2015)

 2

	Variables	AHTS (Brake horsepower: 16,000-20,000 BHP)	PSV (Deck area: 500-900 m²)
		mean	st.dev.	mean	st.dev.
	Dependent variable
	Freight rate (GBP/day)	44,153	46,724	16,584	13,022
	Explanatory variables
	Brake horsepower (BHP)	17,584	1,080
	Deck area (m²)			709	95
	Age (in years)	5.120	4.723	10.347	10.093
	Dynamic positioning class 1 (%)			0.248	0.432
	Dynamic positioning class 2 (%)			0.548	0.498
	Remotely operated vehicle support (%)	0.035	0.185
	Ice class (%)	0.678	0.467	0.054	0.226
	Conventional diesel (%)	0.630	0.483	0.442	0.497
	Activity cargo run (%)	0.146	0.353	0.737	0.440
	rig move (%)	0.658	0.474
	supply (%)			0.248	0.432
	other (%)	0.196	0.397	0.015	0.123
	Contract duration (days)	3.260	3.927	3.057	4.453
	Number of observations	3,109		5,302
	Source: authors' calculations, OSV database. North Sea spot market only.

Table 3 . Estimates of the log of freight rate (2000-2015)

 3 ' calculations, OSV database. Note: estimates from linear regression models, with robust standard errors in parentheses. Significance levels are 1% (***), 5% (**) and 10% (*). The sample is restricted to the North Sea spot market.

	Variables	AHTS (Brake horsepower:	PSV (Deck area: 500-900
		16,000-20,000 BHP)	m²)
		Model 1	Model 2	Model 1	Model 2
	Constant	8.531***	-1.158	8.529***	7.082***
		(0.078)	(1.582)	(0.045)	(0.257)
	Brake horsepower (BHP in log)		1.017***
			(0.162)
	Deck area (m² in log)				0.229***
					(0.034)
	Age (in years /10)		-0.002		0.024
			(0.055)		(0.015)
	Age² (in years /100)		-0.020		-0.021***
			(0.028)		(0.005)
	Dynamic positioning class 1				-0.047***
					(0.014)
	Dynamic positioning class 2				0.026*
					(0.015)
	Remotely operated vehicle support		0.185***
			(0.061)
	Ice Class		0.026		0.060**
			(0.021)		(0.023)
	Conventional Diesel		0.038**		-0.010
			(0.019)		(0.009)
	Activity: cargo run		-0.208***		-0.173***
			(0.035)		(0.045)
	Activity: rig move		0.156***
			(0.027)
	Activity: supply				-0.083*
					(0.044)
	Contract duration (days in log)		0.043*		-0.023**
			(0.026)		(0.011)
	Monthly fixed effects	YES	YES	YES	YES
	Charterer fixed effects	NO	YES	NO	YES
	Number of observations	3,109	3,109	5,302	5,302
	R²	0.730	0.782	0.822	0.856
	Source: authors

Table 4 . Estimates of the distance between Clarksons and transaction-based market indices

 4 ' calculations, OSV database. Note: estimates from linear regression models, with robust standard errors in parentheses. Significance levels are 1% (***), 5% (**) and 10% (*). The sample is restricted to the North Sea spot market.

	Variables	AHTS (16,000-20,000 BHP)		PSV (500-899 m²)
		Model 1	Model 2	Model 3	Model 1	Model 2	Model 3
	Constant	-2,585.493*** -2,562.437*** -2,452.654*** -973.136*** -697.780*** -930.618***
		(453.785)	(447.991)	(405.802)	(109.060)	(91.377)	(104.482)
	Average of freight rate (GBP)	271.473***	279.915***	287.915***	110.732***	91.522***	124.396***
		(44.549)	(46.233)	(43.951)	(11.474)	(9.802)	(12.182)
	St. deviation of net freight rate	80.062	58.393	90.552	46.639	-145.359	-154.427
		(95.934)	(100.481)	(99.317)	(99.549)	(95.893)	(108.671)
	Monthly number of transactions AHTS		-6.233**
			(3.068)
	Monthly number of transactions PSV					-2.029***
						(0.342)
	Monthly use rate of AHTS			-719.837***
				(151.699)
	Monthly use rate of PSV						-244.957***
							(47.089)
	Number of observations	192	192	192	192	192	192
	R²	0.363	0.372	0.445	0.451	0.507	0.496
	Source: authors

The fixture materializes a contract between the seller (vessel owner) and the buyer of the transportation service (charterer) at an agreed price (the freight rate), with the parties usually matched by an intermediary.

See for instance http://splash247.com/bhp-billiton-bypasses-brokers-new-online-auction-platform/.

Without loss of generality, we will normalize the initial value 𝕊𝕊 1 to 1 (or 100) in our empirical analysis.

We also know the identity of the ship owner for each contract but we do not account for this information as most ship owners have very few offshore vessels.

Also, the estimation of (2) on the whole period constrains the various coefficients 𝛽𝛽 to be the same over the whole period. However, the implicit price of some vessel characteristics may change over time. For instance, this is likely to be the case for fuel efficiency in a context where oil price tends to increase.

For details, see https://login.ods-petrodata.com/.

We note that there is some positive correlation between dayrates in the spot/term and AHTS/PSV markets, notably 0.362 between the spot and term AHTS markets, 0.546 between the spot and term PSV markets, 0.776 between the spot AHTS and PSV markets and 0.475 between the term AHTS and PSV markets.

The main charterers for AHTS vessels are Statoil (235 fixtures), BP plc (200 fixtures), Team (173 fixtures) and Nexens Inc. (131 fixtures). The main charterers for PSVs are Team (502 fixtures), BP plc (337 fixtures), Maersk (314 fixtures), ASCO (284 fixtures) and Statoil (275 fixtures).

The main operators are Viking Supply Ships, GulfMark Offshore, Maersk Supply Service, Solstad Shipping and Simon Mokster Shipping for AHTS, Gulf Offshore, Sealion Shipping, GulfMark Offshore, Viking Supply Ships and COH Offshore AS for PSV.

In our dataset, ports where OSVs wait between contracts are Bergen and Stavanger in Norway and Aberdeen in the UK(ICS, 2011). While delivery place is a common variable for deep-sea shipping, this is not the case in the offshore markets, presumably as the services are local (for instance, a voyage from the Norwegian side of the North Sea to the UK side will take less than 24 hours). In the few cases where the loading base is far from the previous port of the vessel, there will be an additional "mobilization" fee included in the contract.

Clarkson Research (2017) also produces these indices at the weekly level, see www.clarksons.net.

Dick, van D., H.E. Haralambides, and A.W. Veenstra (1998): 'Modelling of ocean charter rates; an investigation into non-linearity and parameters '. In: 8th World Conference on Transport Research (WCTR), Antwerp, Belgium, July 1998. Eskandari, H., andE. Mahamoodi (2016): 'A simulation-based multi-objective optimization study of the fleet sizing problem in the offshore industry ', Maritime Economics and Logistics, 18, 436-457. Frances, P. and A.W. Veenstra (1997): 'A co-integration approach to forecasting freight rates in the dry bulk shipping sector ', Transportation Research, 31, 447-458. ICS (2011)


Nadine Malhouroux-Gaffet #

Patrice Pajusco
email: patrice.pajusco@telecom-bretagne.eu

Edgard Haddad

Effect of propagation phenomena on MIMO capacity of wireless systems between 3 and 10 GHz

Multiple Input Multiple Output (MIMO) Antenna systems promise high spectral efficiency [1] over multipath channels. But accurate wideband MIMO channels models are required to optimize these new radio access schemes. This paper focuses on the MIMO channel capacity versus different propagation channel parameters such as frequency or antenna type. This study is based on both experimental and simulated results in a residential environment. The simulations have been computed with a full 3D ray tracing tool while measurements have been performed with a UWB MIMO channel sounder.

I. INTRODUCTION

Increasing the system capacity is a constant goal in the mobile communication research area. The promising MIMO technology is still instigating a large of work especially on propagation study. The potential capacity of a wireless system is indeed directly related to the space-time characteristics of the propagation channel. For the last 10 years, efficient channel sounders have been designed. This equipment requires two antenna arrays to estimate directions of departure and arrival at the transmitter and the receiver respectively. Uniform linear arrays (ULA), uniform planar arrays (UPA) or uniform circular arrays (UCA) are often used. However, all of these arrays have drawbacks which can lead to major biases on the estimate of space-time characteristics [START_REF] Foschini | On Limits of Wireless Communications in a Fading Environment when UsingMultiple Antennas[END_REF]. Thus, a comprehensive space-time characterization of the UWB propagation channel is usually not investigated because a complex measurement setup is required. In the first paragraph, this study presents a simple setup which enables the estimation of the full 3D space-time characterization including frequency dependence.

The second part presents a Ray tracing tool "Matrix" that has been optimized with the measurement carried out with the MIMO UWB experiment presented above. The study of the MIMO capacity behaviour in function of the frequency, the antenna pattern, the numbers of reflection, transmission, and diffraction has been completed with the simulated channel. Simulated data allows analysing all these parameters independently. MIMO channel modelling and the MIMO capacity definition based on the singular value decomposition are presented in the third part.

MIMO capacity variations are analysed in function of several parameters as the antenna pattern, the frequency of simulation, the simulation configuration. A final section presents a comparison between simulated and measured MIMO capacities.

II. MIMO UWB EXPERIMENT

A. Description of the Experimental Set-up

Characterizing the propagation channel is the first step toward the development of accurate propagation models. For this purpose, a new UWB experimentation allowing full 3D space-time characterization has been set up [START_REF] Tesserault | Multi-frequencies characterization of building materials: Angular and polarization analysis[END_REF]. The MIMO UWB channel sounder is based on two uniform circular arrays (UCA). One is located at the transmitter and the other at the receiver. Each UCA is built with 5 UWB antennas orthogonally mounted to each other. These UWB sensors present a pattern frequency dependence shown in Fig. 1 and a picture of the 5-axis antennas is depicted on Fig. 2. The transmission channel was estimated with a vector network analyser (VNA). Two LNAs were used to improve the link budget. The measurement is fully automated with Matlab software. The synoptic of the UWB MIMO channel equipment are depicted in Fig. 3. The MIMO configuration consists of measuring all the links between all the transmitting antennas toward all the receiving antennas for each possible location.

B. Environment and Measurement Campaign

The measurement campaign has been carried out in a typical flat built inside the Orange labs premises. The measurement locations are shown in Fig. 4. Thirteen transmitter positions have been chosen to have different radio link configurations for the characterization of the propagation modelling. The direct path between the transmitter and the receiver ranged from 3 to 8 meters and crossed different architecture elements like plasterboard walls, bearing board walls, doors and some furniture.

III. PROPAGATION CHANNEL SIMULATION

A. Ray Tracing Simulation Parameters

The simulations have been performed with a 3D ray tracing software "MATRIX'' based on geometrical optics (GO) and the uniform geometrical theory of diffraction (UTD).

The ray tracing algorithm has been optimized on the whole UWB bandwidth, taking into account dielectric material coefficients previously determined for a large frequency band (2-12 GHz) [START_REF] Landmann | On the Influence of Data Models on Estimated Angular Distributions in Channel Characterisation[END_REF]. Different combinations of the propagation phenomena are taken into account (Reflection, transmission, diffraction). The simulations are performed in a residential environment described by the properties (dimensions, thickness, permittivity, conductivity) of all architectural elements (walls, floor, ceiling, doors, and windows) and the furniture of each room. The computed ray paths are characterized by their 2×2 complex gain, delay, direction of arrival (DoA), direction of departure (DoD). One example of the Power delay Profile is depicted on Fig. 5. The 100 strongest rays have been taken for the computation of the impulse response.

In order to integrate the antenna pattern in the simulations, each element of the 5-axis antenna has been measured in each azimuth and elevation angles in the two polarizations for frequencies from 2 to 16 GHz.

IV. MIMO CAPACITY ESTIMATION

A. Definition

A MIMO system uses several transmitter antennas (n T ) and several receiver antennas (n R ). The received signal y at the antenna i (i=1,…n R ) is the sum of the transmitted signals .{s 1 ,…,s n } multiplied by the complex gain h i,j of the corresponding radio links and the noise n :

1 T n i i j j i j y h s n y s = + = + = ∑ H n 11 1 1 T R R T n n n n h h h h ⎛ ⎞ ⎜ ⎟ = ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ H … n i
is the equivalent baseband noise whose elements are considered as zero mean circular-symmetric complex additive Gaussian noise (AWGN) samples with a variance of σ 2 . H can be decomposed in singular values:

{ } 1 with b H i i diag λ = = H = UΣV Σ
If the transmitter has no channel state information, the simplest way for power allotment is to distribute the available power equally over the n T antennas. In this case, the capacity is given by: min( , )

2 1 log 1 t r n n i i T C n ρ λ = ⎛ ⎞ = + ⎜ ⎟ ⎝ ⎠ ∑ B.

Computation of the MIMO Capacity

The MIMO-UWB experiment described above corresponds to a virtual 180×180 sensors array. The geometry of the array is circular. This experimentation allows an extraction of subarrays of different sizes, spacing, geometry and numbers of sensors. For instance, four successive antennas can be considered as a small uniform linear array. Hence many configurations can be considered. An average capacity can be estimated by averaging the sub-arrays capacities. This provides a representative estimate of the local channel capacity and this, for different sub-array configurations.

In order to take into account the path loss, the computation of the capacity is realized with a normalized channel power transfer function. The Frobenius norm was used. The normalised MIMO channel matrix is given by :

2 , 1 norm ij i j T R h n n = ∑ H H V. RESULTS ANALYSIS

A. Effect of Propagation Phenomena on the Capacity

A previous investigation based on measurement has been carried out in order to analyse the number of reflections, transmissions and diffractions required to characterize the propagation channel. The result of this study shows wideband parameter stabilization in our environment for a given number of phenomena: 10 transmissions, 6 reflections, 2 diffractions. The channel capacity was computed with these different simulation configurations. The results show a capacity increase of more than 30 bps/Hz between the simulations. In agreement with the power delay profile behaviour (Fig. 5), the capacity value has almost no variation between the configurations with 10 transmissions, 6 reflections and 1 or 2 diffractions.

B. Antenna Pattern Effect on Capacity

Some studies have reported the effect of antenna pattern [START_REF] Pajusco | Novel array structure for space-time characterization of the UWB channel[END_REF] on the estimation of the capacity and suggested to design antennas according to the environment. In an UWB configuration, the antenna pattern depends on the frequency band. Hence, the influence of the antenna pattern on the channel capacity has to be investigated. Ray tracing simulations are very convenient because they allow the decorrelation of the antenna pattern and the propagation channel. So, with a set of rays representing a MIMO link, the estimation of the capacity can be computed with different antenna patterns. The results depicted in figure 7 correspond to the T2-R2 radio link. The simulation is performed at 5 GHz but used different antenna diagrams (fd 103 ,fd 105 ,fd 107 ,fd 109 ) corresponding to the patterns of the UWB antenna at respectively the frequencies:3,5,7,9 GHz. In order to compare only the effect of the antenna pattern, the same normalising factor was applied to all simulations.

A maximal difference of 15 bps/Hz is obtained between the pattern fd 103 and fd 105 . The variation is not linear according to the pattern frequencies. This could be explained by the change of aperture in function of the frequency that leads to a change of integrated rays to the MIMO matrix computation.

C. Effect of the Frequency on Capacity

Experimental analysis of the impact of the frequency has been investigated in [START_REF] Intarapanish | IEEE 2007 Interanaional symposyum on microwave, antenna, propagation and EMC Tecnologies for wireless Communications[END_REF]. These results reported a frequency dependence of the capacity but the influence of the antenna pattern has not been considered. To avoid this effect, ray tracing simulations have been performed at different frequencies and the capacity has been estimated with the same antenna pattern. The results are shown in Fig. 8 for the radio link T2-R2. In a first time, the effect of the frequency has been studied for a single MIMO array (3×3 -spacing=3 cm). If the path loss effect is not taken account, the capacities slightly increase with the frequency. This is explained by the enhancement of the decorrelation of the array sensors. But, if the path loss frequency effect is considered, the capacity strongly decreases (25 bps/Hz) while the frequency increases. Some others capacities simulations have been performed with a sensor spacing set to λ/2 and λ (Fig. 9).These results can be compared to the study of Garcia & Al [START_REF] Intarapanish | IEEE 2007 Interanaional symposyum on microwave, antenna, propagation and EMC Tecnologies for wireless Communications[END_REF]. They kept a constant spacing of one wavelength between sensors and show a capacity decrease while the frequencies increase. Our results are in agreement with them, but the increase of capacity is very low (2-4 bps/Hz) especially for a sensor spacing of λ.

D. Capacity Simulation for Different Locations

Ray tracing simulations have been achieved for each measured location. The capacity has been studied with the path loss effect in order to exactly evaluate the difference of the capacity in a residential flat. The Fig. 10 presents a maximum variation of 20 bps/Hz over the flat. All the MIMO matrix have been normalized to the matrix of the first measurement location.

E. Measured and Simulated MIMO Capacities Comparison

The simulations have been compared to the measurement especially for the transmitter location 2 and 6 which correspond to the better and the worst capacity results in the flat. In the both case, the simulation is slightly undervalued (Fig. 11 &12). That could be explained by a loss of rays or a lack of diffusion modelling in the simulation. The link T2R2 is more affected by this simulation default because at this location, the channel is composed of multipaths stronger than the direct path. For the location T6, the direct path is higher than the multipaths and so the channel could be less affected by the phenomena of diffusion, or omited rays. VI. CONCLUSION An analysis of the modelling parameters of the propagation channel has been realized in order to evaluate its impact on the MIMO capacity. The antenna pattern effect has been evaluated and the frequency dependence of the capacity has been analysed. The simulation results show that the frequency impact on the capacity is mainly due to the path loss in a typical residential flat. The comparison of measurement and simulation reveals an undervaluation of the simulation and shows that the ray-tracing tool needs to be improved. Indeed, a diffusion modelling could give more accuracy and improve the space-time characterisation of the propagation channel.

Fig. 1

 1 Fig. 1 UWB sensor farfield pattern in vertical polarization at a-3, b-5, c-7,d-9 GHz.

Fig. 2 5

 2 Fig. 2 5-axis antenna array

Fig. 4 Fig. 5 :

 45 Fig. 4 Virtual uniform circular array and measurement locations in the residential flat.

Fig. 6

 6 Fig. 6 Mean capacity versus SNR at for a 7×7 array according to the number of propagation (R: Reflection, D: diffraction, 10 Transmissions) phenomena or an isotropic antenna diagram

Fig. 7

 7 Fig. 7 Mean capacity at 5 GHz versus SNR of a 7×7 array (spacing= 2cm) in function of antenna pattern.

Fig. 8

 8 Fig. 8 Mean capacity versus SNR of a 3×3 array according with a 3cm antenna spacing and the same antenna pattern

Fig. 9

 9 Fig.9Mean capacity versus SNR of a 3×3 array according to simulation frequency and sensor spacing for one antenna pattern(5 GHz)

Fig. 10

 10 Fig. 10 Simulated capacity versus SNR of a 3×3 array (spacing= 3cm) 5 GHz.

Fig. 11 :

 11 Fig. 11: Mean capacity versus SNR of a 7×7 array at 3,5,9 GHz for the link T6-R2

Fig. 12 :

 12 Fig. 12: Mean capacity versus SNR of a 7×7 array at 3,5,9 GHz for the link T2-R2