David Gerault
email: dagerault@gmail.com

Marine Minier

Christine Solnon
email: christine.solnon@insa-lyon.fr

Constraint Programming Models for Chosen Key Differential Cryptanalysis

In this paper, we introduce Constraint Programming (CP) models to solve a cryptanalytic problem: the chosen key differential attack against the standard block cipher AES. The problem is solved in two steps: In Step 1, bytes are abstracted by binary values; In Step 2, byte values are searched. We introduce two CP models for Step 1: Model 1 is derived from AES rules in a straightforward way; Model 2 contains new constraints that remove invalid solutions filtered out in Step 2. We also introduce a CP model for Step 2. We evaluate scale-up properties of two classical CP solvers (Gecode and Choco) and a hybrid SAT/CP solver (Chuffed). We show that Model 2 is much more efficient than Model 1, and that Chuffed is faster than Choco which is faster than Gecode on the hardest instances of this problem. Furthermore, we prove that a solution claimed to be optimal in two recent cryptanalysis papers is not optimal by providing a better solution.

Introduction

Cryptography ensures properties such as confidentiality, integrity and signature of communications. Cryptanalysis aims at testing whether these properties are actually guaranteed. Whereas public key cryptography relies on hard problems, symmetric key cryptography relies on simple operations that are iterated many times to speed up encryption/decryption. The most important symmetric key primitives are hash functions and ciphers.

Hash functions guarantee integrity by creating a fixed size fingerprint of messages. Many cryptanalytic results have completely broken the standards MD5, SHA-0 and SHA-1 [START_REF] Wang | How to break md5 and other hash functions[END_REF][START_REF] Wang | Efficient collision search attacks on sha-0[END_REF][START_REF] Karpman | Practical free-start collision attacks on 76-step SHA-1[END_REF] by finding collisions, i.e., messages with a same fingerprint. SAT solvers have been used to find collisions [START_REF] Mironov | Applications of sat solvers to cryptanalysis of hash functions[END_REF][START_REF] De | Inversion attacks on secure hash functions using satsolvers[END_REF][START_REF] Legendre | Encoding hash functions as a sat problem[END_REF] and also against the future hash standard Keccak [START_REF] Morawiecki | A sat-based preimage analysis of reduced keccak hash functions[END_REF].

Ciphers guarantee confidentiality by encoding the original message into a different message, using a key, in such a way that the encoded message can further be decoded into the original one. Stream ciphers encode streams "on the fly", whereas block ciphers split the text in blocks which are encoded separately. Different approaches have been proposed for applying CP to cryptanalysing stream ciphers: [START_REF] Soos | Extending sat solvers to cryptographic problems[END_REF] proposes to solve algebraic systems of equations that link together keys and encoded streams; [START_REF] Mouha | Differential and linear cryptanalysis using mixed-integer linear programming[END_REF] uses mixed integer linear programming to compute bounds for the Enocove-128v2 stream cipher. Since the seminal works of [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF]9], several results appeared on block cipher cryptanalysis [START_REF] Mouha | Differential and linear cryptanalysis using mixed-integer linear programming[END_REF][START_REF] Sun | Extending the applicability of the mixed-integer programming technique in automatic differential cryptanalysis[END_REF][START_REF] Sun | Automatic security evaluation and (related-key) differential characteristic search: Application to simon, present, lblock, DES(L) and other bit-oriented block ciphers[END_REF], mostly based on Mixed-Integer Programming.

Overview of the paper. In this paper, we focus on the cryptanalysis against block ciphers proposed in [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF]9] and described in Section 2. The problem is usually solved in 2 steps: In Step 1, bytes are abstracted by binary values; In Step 2, byte values are searched. In Section 3, we describe a first CP model for Step 1, initially proposed in [START_REF] Minier | Solving a Symmetric Key Cryptographic Problemwith Constraint Programming[END_REF]. This model generates many invalid solutions that are filtered out in Step 2 (as initially proposed in [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF]9]). In Section 4, we introduce new constraints that remove most of these invalid solutions. In Section 5, we briefly describe a CP model for Step 2. In Section 6, we evaluate scale-up properties of two classical CP solvers (Choco and Gecode) and a hybrid CP/SAT solver (Chuffed). We show that the new model for Step 1 is much more efficient than the initial model, and that Chuffed is faster than Choco which is faster than Gecode on the hardest instances of the problem. Furthermore, we prove that a solution claimed to be optimal in [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF]9] is not optimal by providing a better solution. Actually, CP allows us not only to solve cryptanalysis problems more efficiently than the dedicated approaches of [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF]9], but also in a safer way as it is easier to check the correctness of a CP model than the correctness of a dedicated program.

Problem Statement

In this Section, we detail the general structure of the AES (Advanced Encryption Standard) block cipher [START_REF]Advanced Encryption Standard[END_REF]. We then describe what a differential attack is and finally introduce the chosen key differential attack model.

AES block cipher

A block cipher is a function E : {0, 1} n × {0, 1} l → {0, 1} n which, given a binary block X (called plaintext) of length n and a binary key K of length l, outputs a binary ciphered text E(X, K) of length n such that X = E -1 (E(X, K), K).

Most of today's block ciphers have an iterated structure: They apply a round function f r times so that E(X, K) = X r with X 0 = X and X

i+1 = f (X i , K i+1 ) for all i ∈ [0; r -1].
Two famous examples of block ciphers are DES (Data Encryption Standard), which was the encryption standard between 1977 and 2000, and AES [START_REF]Advanced Encryption Standard[END_REF] which is the actual standard since 2001. AES ciphers blocks of length n = 128 bits, where each block is seen as a 4 × 4 matrix of bytes, where a byte is a sequence of 8 bits.

Given a 4×4 matrix of bytes M , we note M The length of keys is l ∈ {128, 192, 256}. The number of rounds depends on the key length: r = 10 (resp. 12 and 14) for l = 128 (resp. 192 and 256). In this Section, we describe AES for l = 128.

The round function f of AES uses an SPN (Substitution-Permutation Network) structure and is described in Fig. 1 for r = 2 rounds. Before the first round, AddRoundKey is applied on the original plaintext X and the initial key

K 0 = K to obtain X 0 = ARK(X, K 0 ). Then, for each round i ∈ [0, r -1]:
-SubBytes, ShiftRows and MixColumns are applied on X i to obtain

Y i = M C(SR(S(X i ))), -KeySchedule is applied on K i to obtain K i+1 = KS(K i ), -AddRoundKey is applied on Y i and K i+1 to obtain X i+1 = ARK(Y i , K i+1 ).
These different operations are described below.

SubBytes S. The S operation, also called S-box, is a non-linear permutation which is applied on each byte of X i separately, i.e., for each j, k ∈

[0, 3], S substitutes X i [j][k] by S(X i [j][k]), according to a lookup table.
ShiftRows SR. SR is a linear mapping that rotates on the left by one byte position (resp. 2 and 3 byte positions) the second row (resp. third and fourth rows) of the current matrix S(X i ), i.e., for each j, k ∈ [0, 3]:

SR(S(X i ))[j][k] = S(X i )[(k + j)%4][k]
MixColumns M C. M C is a linear mapping that multiplies each column of the input matrix SR(S(X i )) by a 4 × 4 fixed matrix chosen for its good properties of diffusion (see [START_REF] Daemen | The Design of Rijndael[END_REF]). In particular, it has the Maximum Distance Separable (MDS) property: For each column j ∈ [0, 3], it ensures:

w(SR(S(X i ))[j]) + w(M C(SR(S(X i )))[j]) ∈ {0, 5, 6, 7, 8}
where w is a function which returns the number of bytes different from 0 8 (we note 0 8 the byte composed of 8 bits equal to 0). AddRoundKey ARK. ARK performs a xor operation (noted ⊕) between Y i and subkey K i+1 , i.e., for each column and row j, k ∈ [0, 3],

ARK(Y i [j][k], K i+1 [j][k]) = Y i [j][k] ⊕ K i+1 [j][k]
KeySchedule KS. The subkey at round 0 is the initial key, i.e., K 0 = K. For each round i ∈ [1, r], the subkey K i is generated from the previous subkey K i-1 by applying the key schedule, i.e., K i = KS(K i-1 ). For keys of length l = 128 bits, each subkey K i is a 4 × 4 byte matrix. KS operates on columns:

-It first computes the first column K i [0] from K i-1 as follows: ∀k ∈ [0; 3], K i [0][k] = K i-1 [0][k] ⊕ S(K i-1 [3][(k + 1)%4])
where S is the SubBytes operation. Moreover, r + 1 predefined constants are added to

K i [0][0].
-For the last 3 columns j ∈ {1, 2, 3}, we have:

K i [j] = K i [j -1] ⊕ K i-1 [j]

Differential Cryptanalysis

Differential cryptanalysis was introduced in 1991 [START_REF] Biham | Differential cryptoanalysis of feal and n-hash[END_REF], and aims at evaluating confidentiality by testing whether it is possible to find the secret key within a reasonable number of trials. The idea is to consider plaintext pairs (X, X ) and to study the propagation of the initial difference between X and X while going through the successive rounds. We note δX i the xor difference between the two plaintexts X i and X i obtained after the ith round of the ciphering of X and X , i.e., δX i = X i ⊕ X i , and we say that Let us keep in mind that the round function f is composed of a set of linear operations (SR, M C, ARK) and a non linear operation (S). The linear operations only move differences to other places. Indeed, for every linear operation l ∈ {SR, M C, ARK}, we have l(A ⊕ B) = l(A) ⊕ l(B). So, we can easily predict how differences are propagated from δX i to δX i+1 by these operations.

δX i [j][k] = X i [j][k] ⊕ X i [j][k]
The non linear operation S has to be studied more carefully. As said before, S operates on each byte X i [j][k] separately. Therefore, we need to study how the S-box propagates differences for a pair (A, B) of bytes. To this aim, we evaluate the probability that the output difference S(A) ⊕ S(B) is equal to β when the input difference A ⊕ B is equal to α, where α and β are bytes. This probability is denoted D α,β and is defined by

D α,β = #{(A, B) ∈ {0, 1} 8 ×{0, 1} 8 | (A⊕B = α) ∧ (S(A)⊕S(B) = β)} 256
For example, let us consider an input difference α = 00000001 and an output difference β = 00100000. For the AES S-box, the transition from 00000001 to 00100000 only occurs for 4 couples of inputs, among the 256 possible couples so that D 00000001,00100000 = 4 256 . For the AES S-box, most of the times the transition probability is equal to 0 256 or 2 256 , and rarely to 4 256 . Note that S is a bijection so that A ⊕ B = 0 8 ⇔ S(A) ⊕ S(B) = 0 8 . As a consequence, D 0 8 ,0 8 = 1. In other words, if there is no difference in the input A ⊕ B, then there is no difference in the output S(A) ⊕ S(B).

Then, for each round i ∈ [0; r], we study difference propagation when the 16 bytes of X i and X i pass through the S-box. For each column j ∈ [0; 3] and each row k ∈ [0; 3], we note

δX i [j][k] = X i [j][k] ⊕ X i [j][k] and δSX i [j][k] = S(X i [j][k]) ⊕ S(X i [j][k]
) the difference for the byte at column j and row k before and after passing the S-box, respectively. The probability of obtaining the output difference

δSX i [j][k] when the input difference is δX i [j][k] is given by D δXi[j][k],δSXi[j][k] .
Hence, the probability of obtaining the output difference δX r = X r ⊕ X r after r rounds given an input difference δX = X ⊕ X is:

p 1 (δX r |δX) = r i=0 3 j=0 3 k=0 D δXi[j][k],δSXi[j][k] (1)
We refer the reader to [START_REF] Biham | Differential cryptoanalysis of feal and n-hash[END_REF] for more details.

A first goal of the attacker is to find the values of δX i for i ∈ {0, . . . r} which maximize the probability p 1 . Once done, the attacker retrieves some partial information on the secret key K from the next subkey K r+1 = KS(K r ). To do so, the attacker has to cipher M chosen plaintext pairs (X, X ) to obtain M ciphered pairs (C, C ). From those pairs (C, C ), the attacker deciphers the last round to partially retrieve δX r according to all possible values of some bits of K r+1 . The correct key will be the one for which the optimal value of δX r (that maximizes p 1 ) appears the most frequently. The number M of plaintext pairs required for the success of the attack may be directly computed from p 1 and is equal to c/p 1 for c a small constant as shown in [START_REF] Biham | Differential cryptoanalysis of feal and n-hash[END_REF].

Chosen Key Differential Cryptanalysis

Today, differential cryptanalysis is public knowledge, so modern block ciphers such as AES have been designed to have proven bounds against differential attacks. However, in 1993, E. Biham proposed a new type of attack called related key attack [START_REF] Biham | New types of cryptoanalytic attacks using related keys (extended abstract)[END_REF] that allows an attacker to inject differences not only between the plaintexts X and X but also between the keys K and K (even if the secret key K stays unknown from the attacker) to try to mount more powerful attacks. The goal of the attack stays the same as previously, i.e., try to find some partial information on the secret key K by testing some bits of the last subkey K r on the veracity of the differential relation which happens with probability p 1 .

We note

K i [j][k] and K i [j][k]
the bytes at column j and row k in the subkeys of K and K at round i, δK i

[j][k] the difference between K i [j][k] and K i [j][k], i.e., δK i [j][k] = K i [j][k] ⊕ K i [j][k], and δSK i [j][k] the difference between S(K i [j][k]) and S(K i [j][k]), i.e., δSK i [j][k] = S(K i [j][k]) ⊕ S(K i [j][k]).
As the only bytes of the subkeys that pass through the S-box are those that are at column j = 3, equation ( 1) is modified by multiplying it by D δKi [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF][k],δSKi [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF][k] , for each round i and each line k, i.e., the goal of the attacker is to find the values of δX i and δK i for i ∈ {0, . . . r} which maximize the probability p 2 defined by equation [START_REF] Biham | Differential cryptoanalysis of feal and n-hash[END_REF].

p 2 (δX r , δK r |δX, δK) = p 1 (δX r |δX) * r i=0 3 k=0 D δKi[3][k],δSKi[3][k]
(2)

Two step solving process for chosen key differential cryptanalysis

Two main papers [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF]9] describe results for the chosen key differential cryptanalysis of the AES and propose algorithms for finding initial pairs of plain texts and keys which maximize the probability p 2 . In both papers, the problem is solved in two steps.

First step: Search of binary solutions. In the first step, each unknown δX i is modeled with a 4 × 4 byte matrix, and a binary variable

∆X i [j][k] is associated with every differential byte δX i [j][k].
These binary variables are equal to 0 if their associated differential bytes are equal to 0 8 , i.e.,

∆X i [j][k] = 0 ⇔ X i [j][k] = X i [j][k] ⇔ δX i [j][k] = 0 8
and they are equal to 1 otherwise. We also associate binary variables

∆K i [j][k] and ∆Y i [j][k] with every differential byte δK i [j][k] = K i [j][k] ⊕ K i [j][k] and δY [j][k] = Y i [j][k] ⊕ Y i [j][k],
respectively. The operations that transform δX into δX r (described in Section 2.1 and Fig. 1), are translated into constraints between these binary variables. In this first step, the goal is to find solutions which satisfy these constraints. Solutions of this first step are called binary solutions. Note that during this first step, the SubBytes operation S is not considered. Indeed, the S-box does not introduce nor remove differences, i.e., for all bytes A and B, (A⊕B = 0 8 ) ⇔ (S(A)⊕S(B) = 0 8 ).

Second step: Search of byte solutions. In the second step, we try to transform binary solutions found in the first step into byte solutions. More precisely, for each binary variable

∆X i [j][k], ∆Y i [j][k] or ∆K i [j][k] set to 1 in the binary solution, we search for a byte value δX i [j][k], δY i [j][k] or δK i [j][k] different from 0 8 so
that the AES transformation rules are satisfied. Note that some binary solutions may not be transformable into byte solutions. These binary solutions are called byte-inconsistent binary solution, whereas binary solutions that can be transformed into byte solutions are called byte-consistent binary solutions. Note also that a byte-consistent binary solution may be transformable into more than one byte solution.

Objective function. The goal is to find a byte solution that maximizes probability p 2 of equation ( 2), while being strictly lower than 1 (i.e., there must be at least one difference between the initial plain texts and keys). It has been shown that a byte solution that maximizes probability p 2 also maximizes the number of factors D α,β of equation [START_REF] Biham | Differential cryptoanalysis of feal and n-hash[END_REF] for which α = β = 0 (because D 0 8 ,0 8 = 1 whereas D α,β ≤ 4 256 if (α, β) = (0 8 , 0 8 )). Therefore, we introduce a variable obj which is equal to the number of ∆X i [j][k] and ∆K i [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF][k] variables of equation ( 2) which are set to 1:

obj = r i=0 3 j=0 3 k=0 ∆X i [j][k] + r i=0 3 k=0 ∆K i [3][k]
We add the constraint obj ≥ 1, to ensure that probability p 2 is strictly lower than 1. To find a byte solution that maximizes p 2 , we first have to find byte-consistent binary solutions that minimize the value of obj. Note that there may exist byteinconsistent binary solutions that have a smaller obj value. However, these binary solutions must be discarded as it is not possible to transform them into byte solutions. Finally, among all byte-consistent binary solutions that minimize obj, we have to search for the one that maximizes the actual probability p 2 associated with its best byte solution.

Existing approaches to solve the problem. In [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF], step 1 is solved with a dedicated Branch & Bound approach. In a preliminary study, we have implemented this approach in C programming language. For r = 3 (resp. r = 4), we found the optimal binary solution in about one hour (resp. 24 hours) on a single core PC. In [9], step 1 is solved by performing a breadth-first traversal of a state-transition graph that has about 2 33.6 nodes for a 128 bit key length. The graph needs 60GB of memory and it is pre-computed in 30 minutes on a 12-core computer for r = 5. Using this graph, binary solutions are found in a few seconds. Both in [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF] and [9], it is claimed that the optimal solution for r = 4 rounds has an objective value obj = 13. We shall see in Section 6 that there exists a better solution.

First CP Model for Step 1

In this section, we describe a CP model for the first step described in Section 2.4. This model was initially introduced in [START_REF] Minier | Solving a Symmetric Key Cryptographic Problemwith Constraint Programming[END_REF].

Variables

Let r be the number of rounds, and let l = 128 be the length of the key. We define the following binary variables (see Fig. 1 for an overview of the bytes associated with these variables):

-For each column and row j, k ∈

[0; 3], ∆X[j][k] is the variable associated with the differential byte δX[j][k] = X[j][k] ⊕ X [j][k].
-For each round i ∈ [0; r] and for each column and row j, k ∈

[0; 3], ∆X i [j][k]
and

∆K i [j][k] are variables associated with differential bytes δX i [j][k] = X i [j][k] ⊕ X i [j][k] and δK i [j][k] = K i [j][k] ⊕ K i [j][k], respectively.
-For each round i ∈ [0; r -1] and for each column and row j, k ∈ [0; 3],

∆Y i [j][k] is the variable associated with the differential byte δY i [j][k] = Y i [j][k] ⊕ Y i [j][k].
All these variables are binary variables, which are set to 0 when the associated differential byte is 0 8 and to 1 otherwise.

Constraints

Constraints correspond to the propagation of differences by operations of the round function f . As said before, the non linear operation S does not imply any constraint as it neither introduces nor removes differences. The linear ARK and KS operations involve xor operations. Therefore, we first define a xor constraint. Then, we define constraints implied by the SR, M C, ARK and KS operations.

xor. Let us consider three differential bytes δA, δB and δC such that δA ⊕ δB = δC. If δA = δB = 0 8 , then δC = 0 8 . If (δA = 0 8 and δB = 0 8 ) or (δA = 0 8 and δB = 0 8 ) then δC = 0 8 . However, if δA = 0 8 and δB = 0 8 , then we cannot know if δC is equal to 0 8 or not. When abstracting differential bytes δA, δB and δC with binary variables ∆A, ∆B and ∆C (which only model the fact that there is a difference or not), we obtain the following definition of the xor constraint:

xor(∆A, ∆B, ∆C) ⇔ ∆A + ∆B + ∆C = 1

AddRoundKey. At the beginning of the ciphering process, ARK performs xor operations on ∆X and ∆K 0 to obtain ∆X 0 , i.e.,

∀(j, k) ∈ [0; 3] 2 , xor(∆X[j][k], ∆K 0 [j][k], ∆X 0 [j][k])
Then, for each round i ∈ [1, r], ARK performs xor operations on ∆Y i-1 and

∆K i to obtain ∆X i , i.e., ∀i ∈ [1, r], ∀(j, k) ∈ [0; 3] 2 , xor(∆Y i-1 [j][k], ∆K i [j][k], ∆X i [j][k])
ShiftRows and MixColumns. For each round i ∈ [0, r-1], SR and MC are applied on ∆X i to obtain ∆Y i , and MC ensures MDS (see Section 2.1), i.e.,

∀i ∈ [0, r -1], ∀j ∈ [0; 3], 3 k=0 ∆X i [(k + j)%4][k] + ∆Y i [j][k] ∈ {0, 5, 6, 7, 8}
KeySchedule. KS is applied at each round i to compute ∆K i from ∆K i-1 , and it is composed of xor operations between some columns of the key, i.e.,

∀i ∈ [1, r], ∀k ∈ [0; 3], xor(∆K i-1 [0][k], ∆K i-1 [3][(k + 1)%4], ∆K i [0][k]) ∀i ∈ [1, r], ∀j ∈ [1; 3], ∀k ∈ [0; 3], xor(∆K i [j -1][k], ∆K i-1 [j][k], ∆K i [j][k])

Objective Variable

Finally, we introduce an integer variable obj, whose domain is [1, l 6 ]5 , and we define obj as the number of differential variables on which a non linear S operation is performed, i.e.,

obj = r i=0 3 j=0 3 k=0 ∆X i [j][k] + r i=0 3 k=0 ∆K i [3][k]

Ordering heuristics

As we want to minimize the number of ∆X i [j][k] and ∆K i [j] [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF] variables set to 1, we add a variable ordering heuristic that first assigns these variables, and a value ordering heuristic that first tries to assign them to 0.

Limitations of the first CP model for Step 1

In [START_REF] Minier | Solving a Symmetric Key Cryptographic Problemwith Constraint Programming[END_REF], we evaluated the CP model described in Section 3 (implemented with Choco 3 [START_REF] Prudhomme | An introduction to choco 3.0: an open source java constraint programming library[END_REF]) on two problems: the optimization problem, the goal of which is to find a binary solution that minimizes the value of obj, and the enumeration problem, the goal of which is to find all binary solutions for a given value of obj (corresponding to the optimal one). These very first experimental results showed us that Choco is able to solve these problems up to r = 5 rounds in a reasonable amount of time. Note that it has been shown in [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF] that it is useless to try to solve these problems for more than 5 rounds because no valid characteristics exist beyond this limit. However, solutions for low values of r are used as a basis to build attacks with larger values of r. For example, [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF] shows how to build an attack for r = 12 and l = 192 by combining 2 solutions with r = 4. Hence, it is very useful to find solutions with lower values of r.

For these two problems, binary solutions are not necessarily byte-consistent. In particular, it may happen that the binary solution of the optimization problem is byte-inconsistent. For instance, for r = 3 rounds, the optimal binary solution has a cost of obj = 3 and there exist 512 binary solutions with this cost. However, none of these solutions are byte-consistent: The optimal byte-consistent binary solution has a cost of obj = 5. When solving the enumeration problem with this cost, we find 21, 504 solutions, among which only 2 are byte-consistent. This means that we spend most of the time at generating useless binary solutions which are discarded in the second step because they are byte-inconsistent. Note that approaches proposed by [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF]9] also suffer from the same problem.

In this section, we introduce new variables and constraints that are added to the first CP model described in Section 3. They are used to infer equality relations between differential bytes, and these relations are used to propagate the MDS property of MixColumns at the byte level. They remove most binary solutions that cannot be transformed into byte solutions, thus speeding up the solution process.

Propagation of MDS at the byte level

For each round i ∈ [0, r -1] and each column j ∈ [0, 3], the MDS property of MixColumns (introduced in Section 2.1) ensures:

3 k=0 w(X i [(k + j)%4][k]) + w(Y i [j][k]) ∈ {0, 5, 6, 7, 8}
At differential byte level, this property still holds:

3 k=0 w(δX i [(k + j)%4][k]) + w(δY i [j][k]) ∈ {0, 5, 6, 7, 8}
In the first model, this property is ensured by the constraint:

3 k=0 ∆X i [(k + j)%4][k] + ∆Y i [j][k] ∈ {0, 5, 6, 7, 8}
However, the MDS property also holds for any xor difference between two different columns in two different rounds of the differential byte model:

∀i 1 , i 2 ∈ [0, r -1], ∀j 1 , j 2 ∈ [0, 3], 3 k=0 w(δX i1 [(k+j 1 )%4][k] ⊕ δX i2 [(k+j 2 )%4][k]) + w(δY i1 [j 1 ][k] ⊕ δY i2 [j 2 ][k]) ∈ {0, 5, 6, 7, 8}
To ensure this property (that removes most byte-inconsistent boolean solutions), we introduce new boolean variables, called equality variables: For each pair of differential bytes δA and δB (in δX i , δY i , and δK i matrices), we introduce the boolean equality variable EQ δA,δB which is equal to 1 if δA = δB, and to 0 otherwise. Using these differential byte equality variables, the MDS property between different columns is ensured by the following constraint:

∀j 1 , j 2 ∈ [0, 3], ∀i 1 , i 2 ∈ [0, r -1], 3 k=0 EQ δXi 1 [(k+j1)%4][k],δXi 2 [(k+j2)%4][k] + EQ δYi 1 [j1][k],δYi 2 [j2][k] ∈ {0, 1, 2, 3, 8}

Constraints on equality variables

In this section, we define constraints that hold on equality variables.

Constraints derived from xor constraints. As pointed out in Section 3.2 when defining the constraint xor(∆A, ∆B, ∆C) (where ∆A, ∆B and ∆C are binary variables associated with differential bytes δA, δB and δC, respectively), if ∆A = ∆B = 1, then we cannot know if ∆C is equal to 0 or 1. However, whenever ∆C = 0 (resp. ∆C = 1), we know for sure that the corresponding byte δC is equal to 0 8 (resp. different from 0 8 ), meaning that the two bytes δA and δB are equal (resp. different), i.e., that EQ δA,δB = 1 (resp. EQ δA,δB = 0). The same reasoning may be done for ∆A and ∆B because (δA ⊕ δB = δC) ⇔ (δB ⊕ δC = δA) ⇔ (δA ⊕ δC = δB). Therefore, we redefine the xor constraint as follows:

xor(∆A, ∆B, ∆C) ⇔ ((∆A + ∆B + ∆C = 1)

∧ (EQ δA,δB = 1 -∆C) ∧ (EQ δA,δC = 1 -∆B) ∧ (EQ δB,δC = 1 -∆A))
Constraints to ensure that equality variables define an equivalence relation. Symmetry is ensured by ∀δA, δB, EQ δA,δB = EQ δ B ,δ A and transitivity by ∀δA, δB, δC, (EQ δA,δB = EQ δB,δC = 1) ⇒ (EQ δ A ,δ C = 1)

Constraints that relate equality variables with binary differential variables. For each pair of differential bytes δA, δB such that the corresponding binary variables are ∆A and ∆B, respectively, we have:

(EQ δA,δB = 1) ⇒ (∆A = ∆B) EQ δA,δB + ∆A + ∆B = 0

Constraints derived from KS

The KeySchedule (described in Section 2.1) mainly performs xor operations: At each round i, the first column K i [0] is obtained by performing a xor between bytes of K i-1 [0] and K i-1 [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF]; for the last three columns j ∈ {1, 2, 3},

K i [j] is obtained by performing a xor between K i-1 [j] and K i [j -1]
. Besides these xor operations, all bytes of K i-1 [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF] pass through the S-box before xoring them with

K i-1 [0] to obtain K i [0]. Therefore, each byte of K i , for each round i ∈ [1, r]
may be expressed as a combination of xor operations between bytes of the initial key K 0 , and bytes obtained by applying the S operation on column 3 of rounds j < i. For example (recall that A ⊕ A = 0 8 and 0 8 ⊕ A = A):

K 2 [1][1] = K 2 [0][1] ⊕ K 1 [1][1] = K 1 [0][1] ⊕ S(K 1 [3][2]) ⊕ K 1 [0][1] ⊕ K 0 [1][1] = S(K 1 [3][2]) ⊕ K 0 [1][1]
When reasoning at the differential byte levels, we have

δK 2 [1][1] = δSK 1 [3][2] ⊕ δK 0 [1][1]
where

δSK 1 [3][2] = S(K 1 [3][2]) ⊕ S(K 1 [3][2]
). As S is a non linear operation, we cannot assume that

δSK 1 [3][2] = S(δK 1 [3][2]). Therefore, δSK 1 [3][2]
is a new differential byte. However, there is a finite number of such new differential bytes: for each round i ∈ [0, r] and each line k ∈ [0, 3], we introduce a new differential byte

δSK i [3][k] = S(K i [3][k]) ⊕ S(K i [3][k])
and a new binary variable ∆SK i [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF][k] which is equal to 0 if δSK i [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF][k] = 0 8 , and to 1 otherwise. Note that

∆SK i [3][k] is a redundant variable which is equal to ∆K i [3][k]. So, we add the constraint ∀i ∈ [1, r], ∀k ∈ [0, 3], ∆K i [3][k] = ∆SK i [3][k]
We introduce this redundant variable because at the byte level this equality no longer holds, i.e.,

δK i [3][k] = A ⇒ δSK i [3][k] = S(A) (because S is a non linear operator such that S(A ⊕ B) = S(A) ⊕ S(B) except when A = B)
, and for the V sets defined below we reason at the byte level.

We propose to exploit the fact that each differential byte of K i is the result of a xor between a finite set of bytes. We first use the KS rules defined in Section 2.1 to build, for each i ∈ [1, r], and j, k ∈ [0, 3], the set V (i, j, k) of all differential bytes (coming either from δK 0 or from the set of new differential bytes δSK i ), such that:

δK i [j][k] = δA∈V (i,j,k) δA For example, V (1, 0, 0) = {δK 0 [0][0], δSK 0 [3][1]}.
Note that these sets are computed before the search and do not depend on the initial values of plaintexts and keys.

For each of these sets, we introduce a set variable which contains the corresponding binary differential variables which are equal to 1:

V 1 (i, j, k) = {∆A | δA ∈ V (i, j, k) ∧ ∆A = 1} For example, if ∆K 0 [1][1] = 1 and ∆SK 1 [3][2] = 0, then V 1 (2, 1, 1) = {∆K 0 [1][1]}.
Whenever two differential key bytes δK i1 [j1][k1] and δK i2 [j2][k2] have the same V 1 sets, then we may infer that δK i1

[j1][k1] = δK i2 [j2][k2]. More precisely, we define the constraint: ∀i 1 , i 2 ∈ [1, r], ∀j 1 , j 2 , k 1 , k 2 ∈ [0, 3], (V 1 (i 1 , j 1 , k 1 ) = V 1 (i 2 , j 2 , k 2 )) ⇒ (EQ δKi 1 [j1][k1],δKi 2 [j2][k2] = 1) Also, if V 1 (i, j, k
) is empty (resp. contains one or two elements), we infer that ∆K i [j][k] is equal to 0 (resp. a variable, or a xor between 2 variables). More precisely, we define the constraints: ∀i ∈

[1, r], ∀j, k ∈ [0, 3], V 1 (i, j, k) = ∅ ⇒ ∆K i [j][k] = 0 V 1 (i, j, k) = {∆A} ⇒ ∆K i [j][k] = 1 ∧ EQ δKi[j][k],δA = 1 V 1 (i, j, k) = {∆A, ∆B} ⇒ xor(∆A, ∆B, ∆K i [j][k])
From a practical point of view, V 1 variables are not modeled with set variables, but with vectors of boolean variables. The dimension of these vectors is equal to the number of possible elements in these sets, i.e., 16 + 4(r + 1) (the 16 bytes of K 0 plus the four bytes that pass through an S-box at each round). Each boolean variable V [p] is equal to 1 if the p th element belongs to V 1 (i.e., if the variable associated with the p th element is equal to 1), and to 0 otherwise. For each of these vectors, we introduce an integer variable which is constrained to be equal to the sum of the variables of the vector.

CP model for Step 2

We have implemented in Choco 3 [START_REF] Prudhomme | An introduction to choco 3.0: an open source java constraint programming library[END_REF] the second step that, given a binary solution, searches for the byte-consistent solution with the highest p 2 value (or prove that there is no byte-consistent solution). The CP model for this second step is rather straightforward and mainly uses table constraints to define relations between the input and the output of the S-box function. The key point is to use a variable ordering that first chooses variables associated with the matrix

∆X i such that j,k ∆X i [j][k] is minimal.
The second step is not a bottleneck and is rather quickly solved by Choco. For example, when l = 128, it is solved in 0.41 (resp. 0.42 and 1.26) seconds, on average, when the number of rounds is r = 3 (resp. r = 4 and r = 5), whereas it is solved in 2.26 seconds when l = 192 and r = 8. Therefore, we have not tried to use other solvers for this step.

Experimental evaluation

In this section, we experimentally compare our two CP models for Step 1: Model 1 refers to the first model introduced in Section 3; Model 2 refers to the first model plus the additional constraints introduced in Section 4. These two models are defined with the MiniZinc modelling language [START_REF] Nethercote | Minizinc: Towards a standard cp modelling language[END_REF]. Model 2 is available at http : //gerault.net/resources/CP _AES.tar.gz.

We compare three solvers on these models: Gecode [START_REF]Gecode: Generic constraint development environment[END_REF] and Choco 4 [START_REF] Prudhomme | An introduction to choco 3.0: an open source java constraint programming library[END_REF], which are classical CP solvers, and Chuffed [START_REF] Chu | Chuffed solver description[END_REF], which is a lazy clause hybrid solver that combines features of finite domain propagation and Boolean satisfiability. All solvers are run on a single core and with default parameters 6 , except option -f for Choco 4 (to break ties of the heuristic described in Section 3.4 with the last conflict heuristic). All runs are limited to one hour of CPU time on a 2.5 -3.5 GHz i7-4710MQ processor with 8 GB of memory.

Table 1 sums up the results for a number of rounds r ∈ {3, 4, 5} 7 . For each round, the objective value obj ranges from the largest value such that Model 1 finds no solution to the smallest value such that there exists a byte-consistent binary solution. Table 1 shows us that Model 2 drastically reduces the number of byte-inconsistent solutions: For example, there are more than 9 * 10 7 byteinconsistent solutions with Model 1 when r=4 and obj=11, whereas there is no solution with Model 2. Hence, Model 2 is much more efficient than Model 1.

For both models, the number of choice points is greater for Gecode than for Choco, and for Choco than for Chuffed. However, choices points are handled faster by Gecode than by Choco (probably because Choco is implemented in Java and Gecode in C++), and faster by Choco than by Chuffed (probably due to lazy clause generation overhead). Therefore, Choco is not faster than Gecode on small instances, and Chuffed is not faster than Choco on small or mediumsize instances. For the hardest instance (r=5; obj=17), Chuffed is nearly twice as fast as Choco, which is nearly twice as fast as Gecode.

All solvers are much faster than the Branch & Bound approach of [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF]: Our C implementation of this approach needs 24 hours to find an optimal binary solution when r = 4. They are also faster and much less memory consuming than the approach of [9], that needs 60GB and 30 minutes on a 12-core computer to pre-compute the graph. For example, for r = 5 and obj = 17, Choco and Chuffed need 400MB and 88MB, respectively.

New results for differential cryptanalysis. We have found two byte-consistent binary solutions with obj = 12 for r = 4 rounds, and we have proven the optimality of these solutions by showing that there does not exist another byte-consistent binary solution with an obj value strictly lower than 12. The optimal byte solution (computed in Step 2) when obj = 12 has a probability p 2 = 2 -79 . The optimal byte solution and its associated binary solution are given in Appendix A. This solution is better than the solution claimed to be optimal in [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF] and [9]: In these papers, authors say that the best byte-consistent binary solution for r = 4 has an obj value equal to 13, and that the optimal byte solution has a probability p 2 = 2 -81 .

Discussion and Conclusion

We have introduced a CP model for solving a problem related to the chosen key differential cryptanalysis of AES with keys of length l = 128. This model follows the classical two step solving process of [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF]9]. In Step 1, we abstract bytes with binary values that indicate whether the byte is equal to 0 8 . In Step 2, we search for non null byte values, for each binary value equal to 1. We have defined new constraints (not used in [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF]9]) which allow us to dramatically reduce the number of binary solutions that cannot be transformed into byte solutions. The idea is to keep track of equalities at the byte level to remove byte-inconsistent solutions at the binary level.

In this paper, we have described models for AES-128, with key length l=128. We have also defined MiniZinc models for AES-192, with l=192. At this time, the best solution we obtained for AES-192 concerns 8 rounds and has obj = 19 active S-boxes. We also plan to extend this work to other families of block ciphers, such as Rijndael [START_REF] Daemen | The Design of Rijndael[END_REF] for which the approach of [9] cannot be used because of its exponential memory complexity.

In our model, we use boolean variables to represent equivalence classes defined by byte equalities: For each pair of bytes, we introduce a boolean variable which is set to 1 if the two bytes are equal, and we explicitly add constraints to ensure symmetry and transitivity of the equality relation. Another possibility would have been to use a graph variable (whose nodes are differential bytes, and edges are byte equality relations), and to post an n-clique global constraint on it, as proposed by Fages [7]: This constraint ensures that the graph is composed of n disjoint cliques, where each clique corresponds to an equivalence class. We have not used this constraint in our model, as it is not available in MiniZinc. We plan to investigate the interest of this constraint using Choco.

Finally, the CP model for Step 2 mainly uses table constraints. Some AES operations operate at the bit level (mostly xor operations), and we plan to improve our model by using bit-vector variables and channeling them with integer variables used to model bytes, as proposed in [START_REF] Michel | Constraint satisfaction over bitvectors[END_REF].

Fig. 1 .

 1 Fig.1. AES ciphering process with r = 2 rounds. Each 4 × 4 array represents a group of 16 bytes. Before the first round, X0 is obtained by applying ARK on the initial text X and the initial key K = K0. Then, for each round i ∈ [1, 2], S, SR and M C are applied on Xi to obtain Yi, KS is applied on Ki-1 to obtain Ki, and ARK is applied on Ki and Yi-1 to obtain Xi. Bytes that pass through the S-box are signaled by an S.

  is a differential byte (for each column and row j, k ∈ [0, 3]).

Table 1 .

 1 Comparison of models and solvers, on the enumeration problem. Each line displays: The number of rounds r, the objective function value obj, the number of byte-consistent binary solutions (S), and the results with models 1 and 2 (number of binary solutions (bin), CPU time in seconds (Time) and number of choice points (CP) for Gecode, Choco 4 and Chuffed). We report '-' when Time is greater than 3600.

					Model 1					Model 2
	r obj S bin Gecode		Choco 4		Chuffed	bin Gecode	Choco 4	Chuffed
			Time CP Time CP	Time CP	Time CP Time CP Time CP
	3 2 0	0 0.0 9E1	0.0 4E1	0.0 5E1 0	0.0 9E1	0.1 4E1	0.0 5E1
	3 3 0 5E2	0.1 2E3	0.4 2E3	0.0 7E2 0	0.0 3E2	0.3 2E2	0.1 2E2
	3 4 0 5E3	1.3 2E4	1.8 1E4	0.2 5E3 0	0.2 9E2	0.5 4E2	0.2 4E2
	3 5 2 2E4	6.0 6E4	5.1 5E4	0.9 2E4 4	0.4 2E3	0.6 1E3	0.6 1E3
	4 8 0	0 0.2 2E4	0.6 1E4	0.3 8E3 0	4.6 1E4	4.9 5E3	6.2 4E3
	4 9 0 2E4	7.1 1E5	5.4 7E4	1.4 4E4 0	8.1 2E4	7.8 8E3	10.7 7E3
	4 10 0 6E6	-	-1161.2 2E7 113.5 6E6 0	14.2 3E4 12.8 1E4	16.2 1E4
	4 11 0 9E7	-	-	-	-1974.5 9E7 0	24.4 5E4 15.5 2E4	25.2 2E4
	4 12 2	-	-	-	-	-	-	-8	44.7 1E5 28.4 5E4	35.7 3E4
	5 10 0	0 1.1 1E5	1.4 5E4	2.3 4E4 0	39.2 3E4 26.8 2E4	37.3 1E4
	5 11 0 3E1 2.0 2E5	2.4 1E5	5.0 7E4 0	63.0 5E4 46.4 3E4	61.5 2E4
	5 12 0 5E5 998.0 4E6 98.3 2E6	48.4 7E5 0 110.0 9E4 74.6 5E4	97.9 3E4
	5 13 0 4E7	-	-	-	-1246.5 5E7 0 187.4 2E5 142.1 9E4 157.6 5E4
	5 14 0	-	-	-	-	-	-	-0 321.7 3E5 247.4 2E5 246.5 8E4
	5 15 0	-	-	-	-	-	-	-10 586.7 5E5 448.2 3E5 408.1 1E5
	5 16 0	-	-	-	-	-	-	-35 1175.8 1E6 770.1 6E5 593.5 2E5
	5 17 6	-	-	-	-	-	-	-50 2879.0 5E6 1524.9 1E6 885.1 4E5

The upper bound l

comes from the fact that D α,β ≤ 2 -6 , ∀(α, β) = (0 8 , 0 8 ), and probability p2 must be larger than 2 -l which corresponds to a probability with uniform distribution of the 2 l possible keys.

We tried other parameter settings. The best results were obtained with default ones.

Let us recall that it has been shown in[START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF] that it is useless to try to solve the problem for more than 5 rounds when the key length is l = 128.

Acknowledgements. Many thanks to Jean-Guillaume Fages, for sending us Choco 4 before the official public release, and to Yves Deville, Pierre Schaus and François-Xavier Standaert for enriching discussions on this work.

A Solution with obj = 12 active S-boxes for AES with r = 4 rounds and l = 128 bits

The byte-consistent binary solution is displayed below. Each binary variable assigned to 1 is colored in blue, and is surrounded in red when it belongs to the objective function (i.e., it passes through an S-box).

The optimal byte solution is displayed below, in hexadecimal notation.

Round δX = X ⊕ X δK = K ⊕ K Init. 0d151846 0dacf2f2 0dfff2f2 0dacf2f2 0 00000000 00ac0000 00000000 00ac0000 0d151846 0d00f2f2 0dfff2f2 0d00f2f2 1 00000000 00ff0000 00ff0000 00000000 0dfff2f2 00ff0000 0dfff2f2 00000000 2 00000000 00ff0000 00000000 00000000 0dfff2f2 0d00f2f2 00000000 00000000 3 00000000 00ff0000 00ff0000 00ff0000 0dfff2f2 00ff0000 00ff0000 00ff0000 End/4 fa000000 faff0000 fa000000 f700f2f2 f7fff2f2 f700f2f2 f7fff2f2 f700f2f2

The corresponding plaintexts X and X and keys K and K are displayed below, in hexadecimal notation. The probability p 2 associated with these plaintexts and keys is p 2 = 2 -79 whereas it is equal to p 2 = 2 -81 in the solution given in [START_REF] Biryukov | Automatic search for related-key differential characteristics in byte-oriented block ciphers: Application to aes, camellia, khazad and others[END_REF] and [9] (solution with obj = 13 active S-boxes).

Round K K 0 00000000 00000000 00000000 00000000 0d151846 0d00f2f2 0dfff2f2 0d00f2f2 1 62636363 62636363 62636363 62636363 6f9c9191 629C6363 6f639191 62636363 2 9b9898c9 f9fbfbaa 9b9898c9 f9fbfbaa 96676a3b f4fb0958 9b9898c9 f9fbfbaa 3 90973450 696ccffa f2f45733 0b0fac99 9d68c6a2 6993cffa f2b5733 0bf0ac99 4 ee60da7b 876a1581 759e42b2 7e91ee2b 19f92889 706ae773 8261b040 89911cd9 Round X X Init. 6b291f8d a800d3d7 f239d5a4 510035ef 663c07cb a5ac2125 ffc62756 5cacc71d 0 6b291f8d a800d3d7 f239d5a4 510035ef 6b291f8d a8acd3d7 f239d5a4 51ac35ef 1 e5000327 00796300 0079005c 0000005a e5000327 00866300 0086005c 0000005a 2 2e2de80b 5186a759 e0d3cbb2 2b02c803 2e2de80b 5179a759 e0d3cbb2 2b02c803 3 5a74f2ae b979ce4a e286aa6a ea86647b 5a74f2ae b986ce4a e279aa6a ea79647b End c501f2fa 4095b6af cdd8f67b 4fadf0a4 3f01f2fa ba6ab6af 37d8f67b b8ad0256


Noëlle Bru
email: nbru@univ-pau.fr

Laurence Despres
email: laurence.despres@ujf-grenoble.fr

Christian Paroissin
email: cparoiss@univ-pau.fr

A comparison of statistical models for short categorical or ordinal time series with applications in ecology

Keywords: time series, categorical variable, ordinal variable, regression model, Markov chain, auto-regressive process, estimation AMS 2000 Classification: 62M10, 62M02, 62M05

We study two statistical models for short-length categorical (or ordinal) time series. The first one is a regression model based on generalized linear model. The second one is a parametrized Markovian model, particularizing the discrete autoregressive model to the case of categorical data. These models are used to analyze two data-sets: annual larch cone production and weekly planktonic abundance.

Introduction

The quest to understand mechanisms behind the temporal dynamics of a natural population (animal or plant) always yields useful information for ecological biodiversity management. The present work was motivated by the analysis of time-fluctuations of an ecological time series, the annual larch cone production. Because of the impracticability of quantitative evaluation of such population size, only semi-quantitative data are often available numerically. Data are coded with finite ordered categories or levels. From this some natural questions arise. Among them, the following one is of crucial interest here: do lagged values determine future production? The same basic problems occur as for classical quantitative time series but here the greatest difficulty stems from the nature of the studied process (as previously mentioned, working on categorical variables induces many difficulties since most of the notions used for quantitative variables have no more sense in such context).

Statistical models have been useful instruments for testing hypothesis concerning the mechanisms behind temporal evolution and to characterize temporal patterns. Two models are used throughout this paper to achieve this goal: a regression model (Fokianos andKedem, 2002, 2003) and a parametrized Markovian model (Jacobs and Lewis, 1978a). The first one is a regression model for categorical time series which is based on generalized linear regression theory [START_REF] Mccullagh | Generalized linear models[END_REF]. Such model extend linear models to accommodate both non-normal response distributions (which is the case in the study of categorical data) and transformations to linearity. So, applying a generalized linear models consists in two choices: a family of probability distribution and a link function between the response and the predictors. For categorical data some widely used models are: multinomial-logit [START_REF] Agresti | Categorical Data Analysis[END_REF] and cumulative odds models [START_REF] Mccullagh | Regression models for ordinal data (with discussion)[END_REF]. Adaptation of such models to categorical times series is easy to do putting past observations at different lags as categorical predictors of the response at time t (Fokianos andKedem, 2002, 2003). The second one is indeed an adaption of the discrete auto-regressive (DAR) model introduced by Jacobs and Lewis (1978a) to the case of categorical time series. As noticed by [START_REF] Mckenzie | Discrete variate time series[END_REF], DAR models ¡¡ would be more suited to modelling dependent sequences of categorical observations, but this does not seem to have been attempted yet ¿¿. To the best of our knowledge, no advance in this direction is made since the paper of McKenzie. These two models have some advantages and some disadvantages which are not necessary the same, implying a complementarity between these two approaches. Among the common advantages, the main one is that they are easy to be interpreted by the practitioners. Since most of the time series in ecology are short-length (for a statistical purpose), we have to consider only models involving a reasonable number of parameters. That is the reason why we will focus on one order lagged model (even these models can be extended easily to large order lag values). Among the inconvenient of the DAR model, the main one is the stationarity of the time series, which can not be checked by any statistical tests (see [START_REF] Mcgee | Coping with nonstationarity in categorical time series[END_REF] for a discussion about several notions of stationarity for categorical time series). However it allows us to derive a simple model for taking into account missing values (a contrario to the regression model, the DAR can not treat directly the case of missing values). Our approach differs highly of the one recently proposed by [START_REF] Bandt | Ordinal time series analysis[END_REF]. Indeed he considers a continuous-state, but non-Gaussian, time series and its analysis relies only on the ordinal property of IR. Moreover his methodology requires a long time series, which is not realistic in many real cases.

Motivation of the present work is the analysis of time series of larch cone production data in spatially disjoint locations in order to determine some temporal patterns of larch cone production dynamic at different locations and to discuss some kind of spatial synchrony. Data are detailed in the first section. Next section 3 is devoted to present two regression models: one for categorical time series and one for ordinal time series. These two regression models have been studied by [START_REF] Fokianos | Regression model for time series analysis[END_REF]2003). In section 4 we adapt the one order discrete auto-regressive model to the context of categorical data (the ordinal characteristic is not taken into account in this model). In particular we develop independence tests and estimators of the various parameters of the model. In section 5, we apply these models to two real data sets: the first one deals with annual larch cone production (over 31 years) and the second one with weekly planktonic abundance (during one year). Last section is devoted to conclusion and discussion.

Motivations

The masting is the intermittent synchronous production of seed crops by a plant population [START_REF] Kelly | Mast seeding in perennial plants: why, how, where?[END_REF]. It often shows an evolved strategy related to others environmental masting patterns such as rainfall, temperatures, . . . Thus variability in seed production according to past values is a good descriptor of environmental changes in climate for example. The information arising from the characterization of temporal patterns on such time series could be used to infer role of environmental parameters and other mechanisms [START_REF] Price | Synchrony and travelling waves of larch bud moth? time series analysis with changing scale[END_REF]. The data accounting cone production were registered for 31 years on four valleys located in the Southern French Alps (in the same area of the Alps called "Brianc ¸onnais"). Here we will consider four different sites selected to be at comparable altitudes (ranging from 1800m to 2200m): Ayes (altitude: 2200 meters), Montgenèvre (altitude: 2200 meters), Névache (altitude: 1800 meters) and Prorel (altitude: 1800 meters). Cone production at a given site was roughly estimated at the beginning of the cone development by counting cones along one meter of branch for at least one hundred randomly selected trees. The intensity of larch cone production at any site was then classified into six classes [START_REF] Roques | The larch cone fly in the french alps[END_REF] from no cones (coded 0) to very heavy crop i.e. more than two hundred cones per tree (coded 4). Annual cone production is considered to be the realization of an ordinal time series with values {0, 0.5, 1, 2, 3, 4} corresponding to a scale classification endowed with a natural ordering. Data are plotted on figure 1.

When studying the dynamic of the larch population on each sampling sites, a first step could be to identify temporal patterns of cone production and then to compare each patterns from one site to others to conclude or not at a spatial synchrony on a "short" regional spatial scale [START_REF] Liebhold | Spatial synchrony in population dynamics[END_REF]. However, the observed series in figure 1 do not exhibit obviously the presence of such patterns. The salient features of the series are: no seasonality, high location to location variability with respect of duration and beginning of intensive larch cone production, presence of missing values, . . . However visual remarks should be considered carefully.

Such time series could appear as too short-length for the statistician who generally needs a lot of information to infer on a phenomenon but the data are collected from 1975 to 2005, which corresponds to an entire career of a biologist!

Regression models for categorical and ordinal time series

The model used here is a generalization of classical regression models to the case of time-dependent categorical observations and was studied by [START_REF] Kauffmann | Regression models for nonstationary categorical time series: asymptotic estimation theory[END_REF] (see also (Fokianos andKedem, 2002, 2003) for a good summary of the main theoretical aspects).

Introduction to generalized linear models for qualitative time series

Assume that the observed series is a particular realization of the stochastic process in discrete time {Y t } which will be described below. Values of Y t are supposed to belong to a finite set E = {1, . . . , k} of k ordered or not categories. Because we are interested in temporal dependence between successive observations, we condition on the observed past. For any positive integer l, let us denote by F t-l the σ-field generated by Y t-1 , Y t-2 , . . . , Y t-l . Let Y t = (Y t,1 , . . . , Y t,k-1 ) ′ where Y t,j equals to 1 if the j-th is observed at time t and 0 otherwise. The analysis of time series based on a generalized linear model suppose that the response variable is influenced by its past values which are viewed as predictors influencing the distribution of Y t by way of a transformation of a linear combination.

IE

[Y t |F t-l ] = h Y ′ t-l β ,
where l is the order of the lag time and Y ′ t-l is the covariate matrix containing the lagged values of the response variable until lag l. In the following we will focus on l ∈ {0, 1, 2} (since we aim at treating short-length time series). The vector β is a vector of time-invariant parameters to be estimate which will reflect the intensity of the dependency between the response and its past. Because the response variable is a categorical time series, we have the following relation:

π t,j,l = IE[Y t,j |F t-l ] = IP(Y t,j = 1|F t-l ),
for every j ∈ {1, . . . , k} and every t{1, . . . , n}, where π t,j,l is a transition probability. Let π t,•,l = (π t,1,l , . . . , π t,k-1,l ). Some adequate regression models for categorical data falls in the family of generalized linear models which links vector of transition probabilities of the response vector to the covariate process through the equation:

π t,•,l := π t,•,l (β) = h(Y ′ t-l β) or h -1 (π t,•,l (β)) = Y ′ t-l β . (1)
In other words, the study of having the response Y t = j at time t is equivalent to carry out a regression on covariates which are the lagged values of the categorical response process. This model is also called a Markov regression model for categorical time series. The function h is called the inverse link function and is related to a link function that describes how the mean depends on the linear predictors. For each response distribution there exists a variety of link functions to connect the mean with the linear predictor. The use of a generalized linear model is the choice of a combination of response distribution and a link function.

On the choice of the link function

The link function should adapted to the type of data [START_REF] Fokianos | Regression theory for categorical time series[END_REF]:

• Nominal data: the most commonly used model for categorical (or nominal) data is the multinomial logit model [START_REF] Agresti | Categorical Data Analysis[END_REF]:

π t,j,l (β) = exp(β ′ j y t-l ) 1 + k-1 q=1 exp(β ′ q y t-l ) , (2)
for any j ∈ {1, . . . , k -1}. This equation also defines log-odds ratios relative to π tm by: log π t,j,l π t,k,l = β ′ j y t-l .

(3)

• Ordinal data: since data are ordinal, its is more convenient to model the cumulative probability function of Y t . For ordered categorical time series a reasonable choice of link function is the logistic distribution one which leads to the proportional odds model [START_REF] Mccullagh | Regression models for ordinal data (with discussion)[END_REF]:

h -1 (x) = 1 1 + exp(-x) . (4)
It follows that the link function is:

log P [Y t ≤ j|F t-l ] P [Y t > j|F t-l ] = Y ′ t-l β .
(5)

Parameters estimation and global adequacy criteria

Since the joint distribution of response and covariates is often not easy to establish, the likelihood methods are not applicable to estimate the vector of regression coefficients β. As we are interested in the estimation of the effects of the covariates on the response, we can use the inference theory based on partial likelihood function. The reader can refer to [START_REF] Fokianos | Regression model for time series analysis[END_REF][START_REF] Viennet | Partial likelihood estimation in categorical time series with stochastic covariates[END_REF] for more details and application. The partial likelihood method leads to non linear equations system. Multinomial models were fitted using the function multinom from library section nnet on R. Proportional-odds logistic regression models were fitted using the function polr from library section MASS on R. The vector of parameters of this model β is estimated using an iterative weighted least squares IWLS [START_REF] Chalmers | Statistical Models in S. Wadsworth & Brooks[END_REF][START_REF] Venables | Modern Applied Statistics with S[END_REF].

The analysis of the global adequacy and goodness of fit of such models to the data is discussed using the Akaike's information criterion (AIC) which also allows to compare several models. The values of this criterion depends on the number of model parameters and penalizes models with large number of parameters. Such consideration is important in the study of short time series where the number of parameters can be rapidly equal to the length of the time series. The chosen model is the one which minimizes the value of AIC among the others.

In this preliminary work, no detailed analysis of the residuals of the models is done. Such analysis is important to assess the goodness of fit between the chosen model and the observed data but was not the priority of this paper.

Discrete auto-regressive model and categorical data

The discrete auto-regressive (DAR) model introduced by Jacobs and Lewis (1978a;[START_REF] Jacobs | Discrete time series generated by mixtures ii: Asymptotic properties[END_REF][START_REF] Jacobs | Discrete time series generated by mixtures iii: autoregressive processes[END_REF] is used here to model categorical data. Some independence tests are developed, using either the Markov property or runs properties. Estimators of the parameters are studied in the precise context of categorical data. Simulated data are used to illustrate numerically the quality of these estimators.

Introduction and model

In a series of papers, Jacobs and Lewis (1978a;[START_REF] Jacobs | Discrete time series generated by mixtures ii: Asymptotic properties[END_REF][START_REF] Jacobs | Discrete time series generated by mixtures iii: autoregressive processes[END_REF] introduced and studied time series models for discrete variables. Among them we will focus here on the discrete auto-regressive of order 1, denoted by DAR(1). Such process {X t } is a discrete-time stochastic process with values on a finite ordered set E = {1, . . . , k} and is defined as follows:

∀t > 0 , X t = V t X t-1 + (1 -V t )Z t ,
where {V t } is a sequence of iid Bernoulli random variables with parameter α ∈ [0; 1] and {Z t } is a sequence of iid random variables having the distribution π on E, the two sequences being independent. Moreover we will assume that X 0 is distributed according to the distribution π, implying that the process {X t } is stationary. The case of α = 1 is not interesting since X t = X 0 , with probability 1, for any t. The case of α = 0 means that the process {X t } is simply a sequence of iid random variables having distribution π. Hence the parameter α could be interpreted as follows: the nearest to 0 α is, the more ¡¡ independent ¿¿ the sequence {X t } is. Indeed, for all h ∈ IN, ρ(h) = α h is the auto-correlation function of a DAR(1) process. Hence DAR(1) models can be used to describe a situation of short range dependency with high correlation. It is easy to prove that stochastic process {X t } as defined above is a Markov chains on E with transition matrix P given by the following equation:

P = αI + (1 -α)Q , where t Q = [ t π| • • • | t π].
Such Markov chain admits obviously a unique stationary probability distribution which is π. One can easily deduce the h-th power of Q and P : for all h

1, Q h = Q and P h = α h I + (1 -α h )Q,
illustrating one more times the role of α.

This stochastic process could be generalized to higher order leading to the DAR(p) model. In fact, these models appear themselves to be a special case of mixture transition distribution (MTD) model introduced by [START_REF] Raftery | A model for high-order markov chains[END_REF]. Thus DAR(p) can be viewed as an alternative to MTD model. According to Raftery, a MTD model fits better data in general than a DAR(p) one, especially for p 3. However here we will prefer to use a DAR(1) model since it has the following advantages over the MTD model: 1) the two parameters α and π play different roles: α is related to the correlation whereas π is the stationary distribution; 2) these models involve generally a reasonable number of parameters (more parsimonious) especially when few data are available; 3) parameters could be easily interpreted by a practitioner. But the special case of DAR(1) model presents the disadvantage of being restrictive over the transition matrix.

Here we are interested on the use of such stochastic processes for modeling categorical variables (here the k different modalities are encoded by using the k first positive integers). It implies that many characteristics of these processes have no sense in such context, as it is the case for the auto-correlation function (see above). Thus estimators developed by Jacobs and Lewis (1983, see pages 28-30) cannot be used. Hence we address the statistical problem of estimating the parameters in a DAR(1) model in presence of categorical data. Assume we observe X 0 , . . . , X n for a fixed value n > 0. First we will test whether {X t } is a sequence of independent random variables (α = 0) or not. In a second step we will estimate all the parameters of the model: α and π. Then we propose a very simple model in order to consider the case of missing observations.

Independence tests

In this section we aim at testing whether {X t } is a sequence of independent random variables (α = 0) or not. Two ways will be investigated. The first one will use the Markov property of the DAR(1) model and the second one will be based on runs property. Anyway, all along this section, the null hypothesis H 0 will be ¡¡ α = 0 ¿¿ and the alternative hypothesis will be ¡¡ α = 0 ¿¿.

χ 2 test based on the Markov property The following test is a classical test for Markov chain (see [START_REF] Reinert | Probabilistic and statistical properties of words: an overview[END_REF] for an illustration in DNA analysis context). We only use the fact that {X t } is a Markov chain, but not the particular structure of its transition matrix. Classical results on Markov chain inference leads to the following estimate for the transition matrix P :

P j,j ′ = N n j,j ′ N n j,•
,

where N n j,j ′ = n i=1 1 1 {Xi-1=j,Xi=j ′ } and N n j,• = j ′ ∈E N n j,j ′ = n i=1 1 1 {Xi-1=j} .
In other words N n j,j ′ is the number of jumps from state j to state j ′ and N n j,• is the number of visits of state j, in the sequence of observations X 0 , . . . , X n .

The null hypothesis can rephrased as follows: P j,j ′ = P j,• P •,j ′ , for any (j, j ′ ) ∈ E 2 . Under H 0 , the maximumlikelihood estimate of P j,j ′ is:

P j,j ′ = P j,• P •,j ′ = N n j,• n -1 N n •,j ′ n -1 , where N n •,j ′ = j∈E N n j,j ′ = n i=1 1 1 {Xi=j ′ } .
Hence one has to consider the following statistics C 2 :

C 2 = j∈E j ′ ∈E [N n j,j ′ -N n j,• N n •,j ′ /(n -1)] 2 N n j,• N n •,j ′ /(n -1)
.

Theorem 4.1 Under the null hypothesis,

C 2 d ----→ n→∞ χ 2 (k-1) 2 .
Some well-known practical restrictions exist in order to be able to apply this test. As example, one can require that P j,j ′ > 5%, for any (j, j ′ ) ∈ E 2 .

Tests based on runs property

Unfortunately we cannot compute the power of the previous test, that is the reason we will now consider a second family of tests. These tests will be based on runs property of the model. Runs in sequence of iid Bernoulli distributions are studied for a very long time: this problem seems to be considered for the first time by Abraham de Moivre in 1756 (problem LXXIV in his book The Doctrine of Chances). For an historical perspective, see the introduction of the Part I of [START_REF] Mood | The distribution theory of runs[END_REF]. Most of the existing papers deal with the case of Bernoulli random variables, but here we are indeed interested in the general discrete case. Few extensions were made in this direction. To the best of our knowledge, [START_REF] Mood | The distribution theory of runs[END_REF] is the first one who studied it.

A run can be defined as follows: it is a consecutive sub-sequence of identical values in a sequence of random numbers. For any j ∈ E, let us denote by R i j,n the number of non-overlapping j-runs of length i in the sequence X 1 , . . . , X n : R i j,n = |{m ; X m-1 = j , , X m = j , , . . . , , X m+i-1 = j , , X m+i = j}| . Let us now define the number R j,n of j-runs and the total number R n of non-overlapping runs: [START_REF] Mood | The distribution theory of runs[END_REF] obtained the limiting distribution of R n after renormalization. Two cases have be distinguished: k = 2 and k > 2.

R j,n = n i=1 R i j,n , and
R n = j∈E R j,n .
Theorem 4.2 (corollary 5 p. 390 and corollary 3 p. 392 in [START_REF] Mood | The distribution theory of runs[END_REF])

1. If k = 2, R n -2nπ 1 π 2 2 nπ 1 π 2 (1 -3π 1 π 2 ) d ----→ n→∞ N (0, 1). 2. If k > 2, 1 √ n R n -n 1 -j∈E π 2 j d ----→ n→∞ N (0, σ 2 π ), with σ 2 π = j∈E π 2 j + 2 j∈E π 3 j - 3( j∈E π 2 j ) 2 .
One can check that in both cases the asymptotic variance is degenerated if and only if there exists j ∈ E such that π j = 1 (and then, for all j ′ = j, π j ′ = 0), then the variance is degenerated. These convergence results could be used to construct an asymptotic test.

An alternative solution could be to consider the longest run in the sequence X 1 , . . . , X n . Indeed there exists many works dealing with the case of either independent trials or Markovian trials. Let us denote by L n the longest length of all runs in the sequence X 1 , . . . , X n . Using the previous notations, we have: [START_REF] Vaggelatou | On the length of the longest run in a multi-state markov chain[END_REF] studied this random variable in the case of multi-state Markovian trials. It requires that {X t } is an irreducible and aperiodic Markov chain on a finite state space E with transition probability matrix P and unique stationary measure π. The Markov chain induced by a DAR(1) model is irreducible and aperiodic if π > 0 (meaning that all the components of π are strictly positive) and α = 1 (see chapter XV of [START_REF] Feller | An introduction to probability theory and its applications: volume I[END_REF]). Let us define the two following quantities:

L n = max{i ; ∃j s.t. R i j,n > 0} .
ρ = max j∈E P jj and π ρ = j∈E : Pjj =ρ π j .
If ρ < 1, then Vaggelatou proved the following asymptotic result (theorem 1 in [START_REF] Vaggelatou | On the length of the longest run in a multi-state markov chain[END_REF]):

IP(L n -[log 1/ρ n] < x) = exp -n(1 -ρ)π ρ ρ [log 1/ρ n]+x-1 + o(1) , (6)
where [•] denotes the integer part and o(1) means that the residual term is ¡¡ small ¿¿ in regard with n. This result extends the classical one obtained many years ago by [START_REF] Goncarov | On the field of combinatory analysis[END_REF] in the case of iid Bernoulli trials.

In both case, L n -[log 1/ρ n] does not have a limit distribution, but only certain sub-sequence; for instance, theorem 2 in [START_REF] Vaggelatou | On the length of the longest run in a multi-state markov chain[END_REF] gives a case where the sub-sequence converges in distribution to the Gumbel distribution. We will use theorem 1 (and not theorem 2) to construct a third (and last) test since we have not enough observations in real situation. Let us denote by ρ 0 and ρ 1 the value of ρ respectively under the null and the alternative hypothesis. Under the null hypothesis, P will be equal to the matrix Q as defined in the introduction: ∀(j, j ′ ) ∈ E 2 , P jj ′ = π j ′ (the transition probability from state j to state j ′ does not depend on j). So we have that ρ 0 = max{π j ; j ∈ E}: ρ 0 < 1 if and only if, for any j ∈ E, π j < 1. Under the alternative hypothesis,

ρ 1 = max{P jj ; j ∈ E} = max{α + (1 -α)π j ; j ∈ E}: ρ 1 < 1 if and only if α < 1 (it is initially assumed)
and for any j ∈ E, π j < 1. Thus we find the same condition in both cases and this assumption is the same as for the previous test. From now we will assume that π > 0 in addition to the previous assumptions (let us recall that we already assume that α = 1). Using theorem 1 of [START_REF] Vaggelatou | On the length of the longest run in a multi-state markov chain[END_REF], one could obtain an asymptotic confidence interval with a prescribed confidence level ε ∈ (0; 1):

IP H0 Ln ∈ Wε,n = log ρ0 - ln(ε/2) n(1 -ρ 0 )π ρ0 ; log ρ0 - ln(1 -ε/2) n(1 -ρ 0 )π ρ0 = 1 -ε ,
Ln = L n -1 (notice that it is corresponding sometimes to the definition of runs: see for instance (Jacobs and Lewis, 1978a)). It follows that the power Π ε of this test is:

Π ε = 1 + IP H1 Ln < log ρ0 - ln(ε/2) n(1 -ρ 0 )π ρ0 -IP H1 Ln < log ρ0 - ln(1 -ε/2) n(1 -ρ 0 )π ρ0 .
To compute Π ε , one has to use equation ( 6) above.

Parameters estimations

The two parameters π and α of such DAR(1) model could be estimated separately since by construction they play different role.

Estimations of π For any j ∈ E and for any i ∈ {1, . . . , n}, let Z ij = 1 1 {Xi=j} : these random variables have the Bernoulli distribution with parameter π j . A natural unbiased estimator of π j is therefore:

π j = 1 n n i=1 Z ij .
Moreover, using the expression of P h given in the introduction, one can easily derive the variance of π j :

Var[ π j ] = 1 n π j (1 -π j ) + 2 n 2 (1 -π j )π j V n (α) ,
and the covariance between π j and π j ′ (with j ′ = j):

Cov[ π j , π j ′ ] = - 2 n 2 π j ′ π j V n (α) , with V n (α) = n h=1 (n -h)α h .
Applying formula (0.113) in (Gradshteyn and Ryznik, 1965) (arithmeticogeometric progression), we obtain the following expression for V n (α):

V n (α) = n -α n 1 -α - α(1 -α n-1 ) (1 -α) 2 -n .
This leads to the following limit for the variances and the covariances:

lim n→∞ nVar[ π j ] = 1 + α 1 -α π j (1 -π j ) ,
and:

lim n→∞ nCov[ π j , π j ′ ] = - 2α 1 -α π j ′ π j .
As a consequence of Bienaymé-Chebychev inequality, the asymptotic result on the variance implies that π j is consistent:

Proposition 4.1 For any α ∈ [0, 1), π j P r ----→ n→∞ π j .
Moreover one can prove the following central limit theorem for π j by application of the ergodic theorem for Markov chain [START_REF] Jones | On the markov chain central limit theorem[END_REF] and Slutsky theorem:

Theorem 4.3 For any α ∈ [0, 1), √ n π j -π j π j (1 -π j ) d ----→ n→∞ N 0, 1 + α 1 -α .
When α = 0, the asymptotic variance equals to 1 as it is well-known for Bernoulli trials. The largest α is, the larger the asymptotic variance is. Since the asymptotic depends on α which is generally unknown, we cannot yet use the last proposition to construct confidence interval. In order to do it, we will need an consistent estimator of α.

Estimations of α We will first consider the maximum likelihood estimator of α, assuming that π is known. Since {X t } is a Markov chain, the log-likelihood is:

L(X 1 , . . . , X n ; α) = (j,j ′ )∈E 2 N j,j ′ log P j,j ′ (α ′ ) ,
where N j,j ′ is defined as in section 2 and where P j,j ′ (α) is the transition probability from state j to state j ′ (which only depends on α). Replacing the expression of transition probabilities, we obtain that:

L(X 1 , . . . , X n ; α) = j∈E   N j,j log(α + (1 -α)π j ) + j ′ ∈E\{j} N j,j ′ log((1 -α)π j ′ )   .
It follows that the maximum likelihood estimator α * 1 of α is the solution of the following equation:

1 n j∈E N n j,j α + (1 -α)π j = 1 .
When π is unknown, one can use the estimation given above and so the plug-in estimator α 1 of α is the solution of the following equation:

1 n j∈E N n j,j α + (1 -α) π j = 1 .
Unfortunately we cannot derive an explicit expression of α 1 . An alternate possible way is to minimize the following function:

Q α = (j,j ′ )∈E 2 ( P j,j ′ -P j,j ′ (α)) 2 .
Solving this optimization problem leads to an explicit expression α 2 * :

α * 2 = j∈E (1 -π j )( P jj -π j ) - j∈E j ′ ∈E\{j} π j ′ ( P jj ′ -π j ′ ) (k -1) j∈E π 2 j + j∈E (1 -π j ) 2 .
It seems to be difficult to establish properties of this intuitive estimator. Hence it is not recommended to use it as an estimator of α. However it could provide a possible initialization for an optimization procedure to obtain a numerical value of α 1 . The corresponding plug-in estimator α 1 of α is given by the following expression:

α 2 = j∈E (1 -π j )( P jj -π j ) - j∈E j ′ ∈E\{j} π j ′ ( P jj ′ -π j ′ ) (k -1) j∈E π 2 j + j∈E (1 -π j ) 2 .
Simulated data The estimators developed above are applied on simulated data in order to evaluate numerically their performance. We simulate data with various values of α, π and n:

• π = ( 1 2 , 1 2 ), π = ( 1 3 , 2 3 ), π = ( 1 3 , 1 3 , 1 3 ) and π = ( 1 4 , 1 2 , 1 4 )
• α ∈ {0.1; 0.2; 0.5; 0.8; 0.9}.

• n ∈ {50; 100; 500}.

Hence two cases are studied: k = 2 and k = 3. For each values of these parameters, we simulate m = 100 independent DAR(1) Markov chain and then we compute the estimators. Unfortunately these we have no guarantee that the two estimators of α belong to the unit interval. Hence we precise for each cases the number of samples on which the computations were done (as it is reasonable, the number of samples is increasing with the number n of observations). Results are given in tables 1 to 4 (only the estimators are computed for simulated data).

A variant with missing observations

Sometimes categorical time series may contain some missing values/observations. Here we now propose a very simple adaptation of the DAR model in order to taking into account the missing values. Since the DAR model is stationary, it will be easy to derive similar expressions as in the initial model.

Assume that at each unit of time, the probability of a missing value equals to β (which does not depend on the time). Hence, if we denote by Z t the values at time t, we have :

Z t = X t w.p. 1 -β -1 w.p. β ,
where -1 is the value corresponding to a missing value and {X t } is DAR(1) stochastic process as described previously. Hence β is the probability that a value is not observed : this probability is assumed to be not depending on t. Since (X t ) is a stationary stochastic process, it follows that (Z t ) is still a Markov chain, but taking values on the set Ẽ = {-1} ∪ E. Its transition probabilities matrix P can be expressed in function of the transition probabilities matrix P of (X t ):

P = 1 -β β t π (1 -β)1 k βP ,
where 1 k is the unit vector of IR k . There is now three parameters to be estimated. Indeed π and α (with the maximum likelihood method) can be estimated as previously. The extra parameter β can be simply estimated as follows:

β = 1 n n i=1 1 1 {Xi=-1} .

Applications to ecological data

Finally tests and estimators are applied in this section to real data. We apply the two models described previously to two real data sets. The first one deals with larch cone production (see section 2) and the second one to planktonic abundance.

Larch cone production

The total number of observations was n = 31, with k = 6 categories. The goal of this work is to study masting on such trees. For a full description of the data, see section 2.

First we apply regression models. Tables 5 and6 contain values of AIC respectively for the categorical time series regression model and the ordinal one. We also indicate the number of parameters to estimate and the number of observations (these time series may contain missing values). For the first case, the independence assumption leads to the better for all sites. For the second case, model with a lag order 1 and 2 fits better for the site Ayes 2200 and Montgenèvre 2200 while the model with only a lag of order 1 fits better for Névache 1800 and the independent model fits better for Prorel 1800. Comparing values of AIC, the model for ordinal time series seems to be more accurate for these data sets.

Second we apply the DAR model. Table 7 contains the estimations of the three parameters for each of the four data sets. In all cases, the independence hypothesis is rejected with the two first tests, while the third one leads to accept the assumption of independence (all with a first type error at 5%). However the power of this last test is more or less weak in all cases (ranging from above 26% to 51%). Thus it is reasonable to reject the assumption of independent observations. In all cases the categorical time series regression model leads to accept the temporal independence between observations. It may be due to the fact that the parameter α in the DAR model is closed to zero. However this parameter can be assumed to be significantly different of zero, according to the performed tests. It is in concordance with the fact that observations are time-dependent when using the ordinal time series regression model. Time series studied here are very short-length and it may the cause that the conclusions based on the ordinal time series regression model and the ones based on the DAR models. Since few data are available, one should prefer to use the DAR models (because it involves less parameters than the other models).

Planktonic abundance

We now consider weekly planktonic (Thalia democratica) abundance data. Data were kindly given by F. Ménard.

In such context, the objective is to test and to compare the temporal patterns from one year to another. Hence we apply the two models described above for four years (1987 to 1990). It follows that each data-set is made of n = 52 observations. Abundances were determined semi-quantitatively according to classes defined on scale of 5

values, E = {1, 2, 3, 4, 5}. The observed series are shown in figure 2 and exhibit the same problems as the larch cone production ones. Notice that the fifth category were not observed for any year. For a complete description of the data, the reader could reefer to [START_REF] Ménard | A stochastic model for ordered categorical time series. application to planktonic abundance data[END_REF] (see also [START_REF] Viennet | Partial likelihood estimation in categorical time series with stochastic covariates[END_REF]).

First we apply regression models. Tables 8 and9 contain the number of parameters to be estimated, the values of AIC and the number of observations, respectively for the categorical time series regression model and the ordinal one. With the model for categorical time series, the model with one order lag fits better all the four years. With the model for categorical time series, model with two order lag fits better for the year 1987 and 1989 while model with only a one order lag fits better for the two other years, 1988 and 1990.

Second we apply the DAR model. Table 10 contains the estimations of the three parameters for each of the four data sets. The two first tests leads to reject the null hypothesis, i.e. to reject the independence of the observations. The third test leads also to reject the independence assumption for the two last year (1989 and 1990) while the null hypothesis is accepted according to this last test for the years 1987 and 1988 (respectively with a power equal to 44% and 66%). Thus it is also reasonabke to reject the assumption of independent observations. We can also analyze these data sets as one unique time series (notice that it was not possible for the previous data-set). Hence the sample size is now n = 208. All the tests lead to reject the assumption of independent observations (with a power equal to 69% for the last one). The last line of table 10 contains the estimations of the three parameters for the whole period.

Here the situation is totally different than previously. Indeed the categorical time series regression model and the DAR model lead to the same conclusion, i.e. a one order lag dependence in the time series. One can notice that for these data sets the parameter α (of the DAR model) is now between 0.308 and 0.540. However the ordinal time series regression model leads in some case to a two order lag dependence. Since the number of observations is almost the twice than for the first data sets, one should rather prefer to use the ordinal time series regression model.

Conclusions and discussions

Applications to real data achieve to convince that these two complementary models are relevant for practical purpose. Based on the result obtained over real data, one can conclude that either the ordinal time-series regression model or the DAR model should be used to treat such data. Indeed in all cases the categorical time-series regression model seems not to present advantages over the two other models. The choice between the ordinal time-series regression model and the DAR model depends highly on the context, i.e. essentially on the number of observations and the number of parameters to be estimated. [START_REF] Fokianos | Regression theory for categorical time series[END_REF] claimed that ¡¡ the regression methodology can discover dependencies in the DNA sequence data which cannot assessed by a Markov model ¿¿. However data treated here can serve as a counter-example of this sentence. For the two data sets studied here, conclusions based on the regression models and the one based on the Markov model are almost identical.

From this work, one can conclude that in any case the DAR model has only few parameters to be estimated, but with an equal to number of unknown parameters (as in the example of larch cone production) one has to prefer the ordinal time-series regression model.

However both suffers of relying on assumptions or simplifications. Hence these models could be extended in the following ways:

• Stationarity: the DAR model is strongly stationary (in the sense defined by [START_REF] Mcgee | Coping with nonstationarity in categorical time series[END_REF]).

This assumption should be checked with any statistical tests. However no test of stationarity of a categorical time series has been developed to the best of our knowledge. Anyway when dealing with short-length time series stationary assumption is not really restrictive. Otherwise a solution could be in applying the de-trend algorithm suggested by [START_REF] Mcgee | Coping with nonstationarity in categorical time series[END_REF]. However their algorithm is more and more computationally complex as the number of states is increasing (in fact they mainly consider the binary case).

We do not focus here on the study of the possible stationarity of a categorical time series which will be done in a future work. A major advantage of regression model is that it is not necessary to have stationarity.

• Higher dependence order and number of parameters: since we consider the case of short-length data, we limit our study to one or two order lagged models. Indeed both models could be applied to p-th order lagged models. However, for the regression model, a large value of p implies a large number of parameters to estimate, that may induce some numerical instability (due to correlation between the regressors). The number of parameters in a DAR(p) model is lower, but if p > 1 we have no more the Markov property.

• Environmental factors: these models do not include environmental covariates. The regression model could easily integrate such situations, as shown by [START_REF] Fokianos | Regression model for time series analysis[END_REF]2003). For the DAR model, it is not so easy. A solution could be to consider inhomogeneous Markov chain or a state-space model.

Biographical sketches

Figure 1 :

 1 Figure 1: Annual larch production in four sites in the Southern Alps

  Noëlle Bru is an Assistant Professor at the IUT STID of Université de Pau et des Pays de l'Adour and a Researcher at the Laboratoire de Mathématiques Appliqués de Pau. Since her PhD thesis, her topics of interest is both theoretical and applied statistics with emphasis to environmental data analysis.Laurence Despres is an Assistant Professor at the Université Joseph Fourier and a Researcher at the Laboratoire d'Ecologie Alpine. Her main research interests are in the evolutionary ecology of species interactions and coevolution.Christian Paroissin is an Assistant Professor at the Département Sciences et Techniques of Université de Pau et des Pays de l'Adour and a Researcher at the Laboratoire de Mathématiques Appliqués de Pau. His topics of interest is applied probability and statistics with interest to applied contexts (engineering, theoretical computer science, biology, . . . ).

Table 1 :

 1 Results obtained with π = ( 1

	α	n	π	α 1	m 1	α 2	m 2
		50	(0.509;0.491) 0.140 65	0.149 72
	0.1	100 (0.499;0.501) 0.114 82	0.121 85
		500 (0.503;0.497) 0.092 100 0.094 100
		50	(0.511;0.489) 0.174 83	0.183 90
	0.2	100 (0.510;0.490) 0.180 97	0.192 97
		500 (0.499;0.501) 0.191 100 0.194 100
		50	(0.501;0.499) 0.451 99	0.481 99
	0.5	100 (0.510;0.490) 0.465 100 0.479 100
		500 (0.499;0.501) 0.494 100 0.497 100
		50	(0.543;0.457) 0.711 99	0.753 99
	0.8	100 (0.534;0.466) 0.765 100 0.784 100
		500 (0.502;0.498) 0.795 100 0.799 100
		50	(0.585;0.415) 0.754 99	0.798 99
	0.9	100 (0.539;0.461) 0.857 100 0.882 100
		500 (0.503;0.494) 0.890 100 0.894 100
						2 , 1 2 )

Table 3 :

 3 Results obtained with π = (1

	3 , 1 3 , 1 3 )

Table 4 :

 4 Results obtained with π = ( 1

	4 , 1 2 , 1 4 )

Table 5 :

 5 Categorical time-series regression models applied to annual larch cones production

Acknowledgments

We wish to thank Patricia Jacobs (Naval Postgraduate School, Monterey, California, USA), Ian R. Harris (Southern Methodist University, Texas, USA), Alain Latour (LabSAD, Grenoble, France), Monnie McGee (Southern Methodist University, Texas, USA) and Frédéric Ménard (IRD, Montpellier, France). Thanks also to the numerous scientists and crew members who conducted the experimental campaigns because our analysis is based on their hard work. This work is a contribution to the understanding of larch cone production into an IFB (Institut Franc ¸ais de la Biodiversité) project.