Daniil Ryabko
email: daniil@ryabko.net

Characterizing predictable classes of processes

The problem is sequence prediction in the following setting. A sequence x1, . . . , xn, . . . of discrete-valued observations is generated according to some unknown probabilistic law (measure) µ. After observing each outcome, it is required to give the conditional probabilities of the next observation. The measure µ belongs to an arbitrary class C of stochastic processes. We are interested in predictors ρ whose conditional probabilities converge to the "true" µ-conditional probabilities if any µ ∈ C is chosen to generate the data. We show that if such a predictor exists, then a predictor can also be obtained as a convex combination of a countably many elements of C. In other words, it can be obtained as a Bayesian predictor whose prior is concentrated on a countable set. This result is established for two very different measures of performance of prediction, one of which is very strong, namely, total variation, and the other is very weak, namely, prediction in expected average Kullback-Leibler divergence.

Introduction

Given a sequence x 1 , . . . , x n of observations x i ∈ X , where X is a finite set, we want to predict what are the probabilities of observing x n+1 = x for each x ∈ X , before x n+1 is revealed, after which the process continues. It is assumed that the sequence is generated by some unknown stochastic process µ, a probability measure on the set of one-way infinite sequences X ∞ . The goal is to have a predictor whose predicted probabilities converge (in a certain sense) to the correct ones (that is, to µ-conditional probabilities). In general this goal is impossible to achieve if nothing is known about the measure µ generating the sequence. In other words, one cannot have a predictor whose error goes to zero for any measure µ. The problem becomes tractable if we assume that the measure µ generating the data belongs to some known class C. The questions addressed in this work are a part of the following general problem: given an arbitrary set C of measures, how can we find a predictor that performs well when the data is generated by any µ ∈ C, and whether it is possible to find such a predictor at all. An example of a generic property of a class C that allows for construction of a predictor, is that C is countable. Clearly, this condition is very strong. An example, important from the applications point of view, of a class C of measures for which predictors are known, is the class of all stationary measures. The general question, however, is very far from being answered.

The contribution of this work to solving this question is in that we we provide a specific form in which to look for a solution to the general problem. More precisely, we show that if a predictor exists, then a predictor can also be obtained as a weighted sum of a countably many elements of C. This result can also be viewed as a justification of the Bayesian approach to sequence prediction: if there exists a predictor which predicts well every measure in the class, then there exists a Bayesian predictor (with a rather simple prior) that has this property too. In this respect it is important to note that the result obtained about such a Bayesian predictor is pointwise (holds for every µ in C), and stretches far beyond the set its prior is concentrated on.

The motivation for studying predictors for arbitrary classes C of processes is two-fold. First of all, prediction is a basic ingredient for constructing intelligent systems. Indeed, in order to be able to find optimal behaviour in an unknown environment, an intelligent agent must be able, at the very least, to predict how the environment is going to behave (or, to be more precise, how relevant parts of the environment are going to behave). Since the response of the environment may in general depend on the actions of the agent, this response is necessarily non-stationary for explorative agents. Therefore, one cannot readily use prediction methods developed for stationary environments, but rather has to find predictors for the classes of processes that can appear as a possible response of the environment.

Apart from this, the problem of prediction itself has numerous applications in such diverse fields as data compression, market analysis, bioninformatics, and many others. It seems clear that prediction methods constructed for one application cannot be expected to be optimal when applied to another. Therefore, an important question is how to develop specific prediction algorithms for each of the domains. In order to do this, the first step is to understand for which classes of problems (i.e. sets of measures generating the data) a predictor exists.

Prior work. As it was mentioned, if the class C of measures is countable (that is, if C can be represented as C := {µ k : k ∈ N}), then there exists a predictor which performs well for any µ ∈ C. Such a predictor can be obtained as a Bayesian mixture ρ S := k∈N w k µ k , where w k are summable positive real weights, and it has very strong predictive properties; in particular, ρ S predicts every µ ∈ C in total variation distance, as follows from the result of [START_REF] Blackwell | Merging of opinions with increasing information[END_REF]. Total variation distance measures the difference in (predicted and true) conditional probabilities of all future events, that is, not only the probabilities of the next observations, but also of observations that are arbitrary far off in the future (see formal definitions below). In the context of sequence prediction the measure ρ S was first studied by [Solomonoff, 1978]. Since then, the idea of taking a convex combination of a finite or countable class of measures (or predictors) to obtain a predictor permeates most of the research on sequential prediction (see, for example, [START_REF] Cesa-Bianchi | Prediction, Learning, and Games[END_REF]) and some related topics in AI [Hutter, 2005, Ryabko andHutter, 2008a]. In practice it is clear that, on the one hand, countable models are not sufficient, since already the class µ p , p ∈ [0, 1] of Bernoulli i.i.d. processes, where p is the probability of 0, is not countable. On the other hand, prediction in total variation can be too strong to require; predicting probabilities of the next observation may be sufficient, maybe even not on every step but in the Cesaro sense. A key observation here is that a predictor ρ S = w k µ k may be a good predictor not only when the data is generated by one of the processes µ k , k ∈ N, but when it comes from a much larger class. Let us consider this point in more detail. Fix for simplicity X = {0, 1}. The Laplace predictor λ(

x n+1 = 0|x 1 , . . . , x n ) = #{i≤n:xi=0}+1 n+|X |
predicts any Bernoulli i.i.d. process: although convergence in total variation distance of conditional probabilities does not hold, predicted probabilities of the next outcome converge to the correct ones. Moreover, generalizing the Laplace predictor, a predictor λ k can be constructed for the class M k of all k-order Markov measures, for any given k. As was found by [Ryabko, 1988], the combination ρ R := w k λ k is a good predictor not only for the the set ∪ k∈N M k of all finite-memory processes, but also for any measure µ coming from a much larger class: that of all stationary measures on X ∞ . Here prediction is possible only in the Cesaro sense (more precisely, ρ R predicts every stationary process in expected time-average Kullback-Leibler divergence, see definitions below). The Laplace predictor itself can be obtained as a Bayes mixture over all Bernoulli i.i.d. measures with uniform prior on the parameter p (the probability of 0). However, as was observed in [Hutter, 2007] (and as is easy to see), the same (asymptotic) predictive properties are possessed by a Bayes mixture with a countably supported prior which is dense in [0, 1] (e.g. taking ρ := w k δ k where δ k , k ∈ N ranges over all Bernoulli i.i.d. measures with rational probability of 0). For a given k, the set of k-order Markov processes is parametrized by finitely many [0, 1]-valued parameters. Taking a dense subset of the values of these parameters, and a mixture of the corresponding measures, results in a predictor for the class of k-order Markov processes. Mixing over these (for all k ∈ N) yields, as in [Ryabko, 1988], a predictor for the class of all stationary processes. Thus, for the mentioned classes of processes, a predictor can be obtained as a Bayes mixture of countably many measures in the class. An additional reason why this kind of analysis is interesting is because of the difficulties arising in trying to construct Bayesian predictors for classes of processes that can not be easily parametrized. Indeed, a natural way to obtain a predictor for a class C of stochastic processes is to take a Bayesian mixture of the class. To do this, one needs to define the structure of a probability space on C. If the class C is well parametrized, as is the case with the set of all Bernoulli i.i.d. process, then one can integrate with respect to the parametrization. In general, when the problem lacks a natural parametrization, although one can define the structure of the probability space on the set of (all) stochastic processes in many different ways, the results one can obtain will then be with probability 1 with respect to the prior distribution (see, for example, [START_REF] Jackson | Bayesian representation of stochastic processes under learning: de Finetti revisited[END_REF]), while pointwise consistency cannot be assured (see e.g. [START_REF] Diaconis | On the consistency of Bayes estimates[END_REF]). Re-sults with prior probability 1 can be hard to interpret if one is not sure that the structure of the probability space defined on the set C is indeed a natural one for the problem at hand (whereas if one does have a natural parametrization, then usually results for every value of the parameter can be obtained, as in the case with Bernoulli i.i.d. processes mentioned above). The results of the present work show that when a predictor exists it can indeed be given as a Bayesian predictor, which predicts every (and not almost every) measure in the class, while its support is only countable.

The results. Here we show that if there is a predictor that performs well for every measure coming from a class C of processes, then a predictor can also be obtained as a convex combination k∈N w k µ k for some µ k ∈ C and some w k > 0, k ∈ N. This holds if the prediction quality is measured by either total variation distance, or expected average KL divergence: one measure of performance that is very strong, the other rather weak. The analysis for the total variation case relies on the fact that if ρ predicts µ in total variation distance, then µ is absolutely continuous with respect to ρ, so that ρ(x 1..n )/µ(x 1..n ) converges to a positive number with µ-probability 1 and with a positive ρ-probability. However, if we settle for a weaker measure of performance, such as expected average KL divergence, measures µ ∈ C are typically singular with respect to a predictor ρ. Nevertheless, since ρ predicts µ we can show that ρ(x 1..n )/µ(x 1..n ) decreases subexponentially with n (with hight probability), and then we can use this ratio as an analogue of the density for each time step n, and find a convex combination of countably many measures from C that has desired predictive properties for each n. Combining these predictors for all n then results in a predictor that predicts every µ ∈ C in average KL divergence. The proof techniques developed have a potential to be used in solving other questions concerning sequence prediction, in particular, the general question of how to find a predictor for an arbitrary class C of measures.

Preliminaries

Let X be a finite set. The notation x 1..n is used for x 1 , . . . , x n . We consider stochastic processes (probability measures) on (X ∞ , F ) where F is the sigmafield generated by the cylinder sets [x 1..n ], x i ∈ X , n ∈ N, where [x 1..n ] is the set of all infinite sequences that start with x 1..n . For a finite set A denote |A| its cardinality. We use E µ for expectation with respect to a measure µ.

Next we introduce the measures of the quality of prediction used in this paper. For two measures µ and ρ we are interested in how different the µ-and ρ-conditional probabilities are, given a data sample x 1..n . Introduce the total variation distance

v(µ, ρ, x 1..n ) := sup A∈F |ρ(A|x 1..n ) -µ(A|x 1..n )|.
Definition 1. We say that ρ predicts µ in total variation if

v(µ, ρ, x 1..n ) → 0 µ-a.s.
This convergence is rather strong. In particular, it means that ρ-conditional probabilities of arbitrary far-off events converge to µ-conditional probabilities. Moreover, ρ predicts µ in total variation if [START_REF] Blackwell | Merging of opinions with increasing information[END_REF] and only if [START_REF] Kalai | Weak and strong merging of opinions[END_REF] µ is absolutely continuous with respect to ρ.

Thus, for a class C of measures there is a predictor ρ that predicts every µ ∈ C in total variation if and only if every µ ∈ C has a density with respect to ρ. Although such sets of processes are rather large, they do not include even such basic examples as the set of all Bernoulli i.i.d. processes. That is, there is no ρ that would predict in total variation every Bernoulli i.i.d. process measure δ p , p ∈ [0, 1], where p is the probability of 0. Therefore, perhaps for many (if not most) practical applications this measure of the quality of prediction is too strong, and one is interested in weaker measures of performance.

For two measures µ and ρ introduce the expected cumulative Kullback-Leibler divergence (KL divergence) as

d n (µ, ρ) := E µ n t=1 a∈X µ(x t = a|x 1..t-1 ) log µ(x t = a|x 1..t-1 ) ρ(x t = a|x 1..t-1 ) , (1)
In words, we take the expected (over data) average (over time) KL divergence between µ-and ρ-conditional (on the past data) probability distributions of the next outcome.

Definition 2. We say that ρ predicts µ in expected average KL divergence if

1 n d n (µ, ρ) → 0.
This measure of performance is much weaker, in the sense that it requires good predictions only one step ahead, and not on every step but only on average; also the convergence is not with probability 1 but in expectation. With prediction quality so measured, predictors exist for relatively large classes of measures; most notably, [Ryabko, 1988] provides a predictor which predicts every stationary process in expected average KL divergence. A simple but useful identity that we will need (in the context of sequence prediction introduced also in [Ryabko, 1988]) is the following

d n (µ, ρ) = - x1..n∈X n µ(x 1..n ) log ρ(x 1..n ) µ(x 1..n ) , (2)
where on the right-hand side we have simply the KL divergence between measures µ and ρ restricted to the first n observations. Thus, the results of this work will be established with respect to two very different measures of prediction quality, one of which is very strong and the other rather weak. This suggests that the facts established reflect some fundamental properties of the problem of prediction, rather than those pertinent to particular measures of performance. On the other hand, it remains open to extend the results below to different measures of performance.

Main results

Theorem 1. Let C be a set of probability measures on X ∞ . If there is a measure ρ such that ρ predicts every µ ∈ C in total variation, then there is a sequence µ k ∈ C, k ∈ N such that the measure ν := k∈N w k µ k predicts every µ ∈ C in total variation, where w k are any positive weights that sum to 1. This relatively simple fact can be proven in different ways, relying on the equivalence of the statements "ρ predicts µ in total variation distance" and "µ is absolutely continuous with respect to ρ." The proof presented below uses techniques that can be then generalized to the case of prediction in expected average KL-divergence, where in all interesting cases all measures µ ∈ C are singular with respect to any predictor that predicts all of them. The idea of the proof of Theorem 1 is as follows. For each measure µ ∈ C we find the set T µ of sequences x 1 , x 2 , . . . on which the density of µ with respect to ρ exists and is non-zero. Such a set has µ-probability 1, and, by absolute continuity, a positive ρ-probability. The idea is then to cover the union ∪ µ∈C T µ with countably many of these sets, and then construct a new predictor as a sum of the corresponding measures. To find this countable collection of sets T µ , we first find a largest (up to an ε 1 ) one with respect ρ, then the one who has a largest (up to an ε 2 ) part not covered by the first set, and so on (where ε k are decreasing). Then we show that any strictly convex combination of the resulting sequence of measures has the property that any measure in C is absolutely continuous with respect to it.

Proof. We break the (relatively easy) proof of this theorem into 3 steps, which will make the (more involved) proof of the next theorem more understandable.

Step 1: densities. For any µ ∈ C, since ρ predicts µ in total variation, µ has a density (Radon-Nikodym derivative) f µ with respect to ρ. Thus, for the set T µ of all sequences x 1 , x 2 , ... ∈ X ∞ on which f µ (x 1,2,... ) > 0 (the limit lim n→∞ ρ(x1..n) µ(x1..n) exists and is finite and positive) we have µ(T µ ) = 1 and ρ(T µ ) > 0. Next we will construct a sequence of measures µ k ∈ C, k ∈ N such that the union of the sets T µ k has probability 1 with respect to every µ ∈ C, and will show that this is a sequence of measures whose existence is asserted in the theorem statement.

Step 2: a countable cover and the resulting predictor. Let ε k := 2 -k and let m 1 := sup µ∈C ρ(T µ ). Clearly, m 1 > 0. Find any µ 1 ∈ C such that ρ(T µ1 ) ≥ m 1 -ε 1 , and let T 1 = T µ1 . For k > 1 define m k := sup µ∈C ρ(T µ \T k-1 ). If m k = 0 then define T k := T k-1 , otherwise find any µ k such that ρ(T µ k \T k-1 ) ≥ m k -ε k , and let T k := T k-1 ∪ T µ k . Define the predictor ν as ν := k∈N w k µ k .

Step 3: ν predicts every µ ∈ C. Since the sets T 1 , T 2 \T 1 , . . . , T k \T k-1 , . . . are disjoint, we must have ρ(T k \T k-1 ) → 0, so that m k → 0. Let

T := ∪ k∈N T k .
Fix any µ ∈ C. Suppose that µ(T µ \T ) > 0. Since µ is absolutely continuous with respect to ρ, we must have δ := ρ(T µ \T ) > 0. Then for every k > 1 we have

m k = sup µ ′ ∈C ρ(T µ ′ \T k-1 ) ≥ ρ(T µ \T k-1 ) ≥ δ > 0,
which contradicts m k → 0. Thus, we have shown that

µ(T ∩ T µ ) = 1. (3)
Let us show that every µ ∈ C is absolutely continuous with respect to ν. Indeed, fix any µ ∈ C and suppose µ(A) > 0 for some A ∈ F. Then from (3) we have µ(A ∩ T ) > 0, and, by absolute continuity of µ with respect to ρ, also ρ(A ∩ T ) > 0. Since T = ∪ k∈N T k we must have ρ(A ∩ T k ) > 0 for some k ∈ N. Since on the set T k the measure µ k has non-zero density f µ k with respect to ρ, we must have µ

k (A ∩ T k ) > 0. (Indeed, µ k (A ∩ T k ) = A∩T k f µ k dρ > 0.) Hence, ν(A ∩ T k ) ≥ w k µ k (A ∩ T k ) > 0,
so that ν(A) > 0. Thus, µ is absolutely continuous with respect to ν, and so ν predicts µ in total variation distance.

Theorem 2. Let C be a set of probability measures on X ∞ . If there is a measure ρ such that ρ predicts every µ ∈ C in expected average KL divergence, then there is a sequence µ k ∈ C, k ∈ N such that the measure ν := k∈N w k µ k predicts every µ ∈ C in expected average KL divergence, where w k are some positive weights.

A difference worth noting with respect to the formulation of Theorem 1 (apart from a different measure of divergence) is in that in the latter the weights w k can be chosen arbitrarily, while in Theorem 2 they can not. In general, the statement " k∈N w k ν k predicts µ in expected average KL divergence for some choice of w k , k ∈ N" does not imply " k∈N w ′ k ν k predicts µ in expected average KL divergence for every summable sequence of positive w ′ k , k ∈ N," while the implication trivially holds true if the expected average KL divergence is replaced by the total variation. An interesting related question (which is beyond the scope of this paper) is how to chose the weights to optimize the behaviour of a predictor before asymptotic.

The idea of the proof is as follows. For every µ and every n we consider the sets T n µ of those x 1..n on which µ is greater than ρ. These sets have to have (from some n on) a high probability with respect to µ. Then since ρ predicts µ in expected average KL divergence, the ρ-probability of these sets cannot decrease exponentially fast (that is, it has to be quite large). (The sequences µ(x 1..n )/ρ(x 1..n ), n ∈ N will play the role of densities of the proof of Theorem 1, and the sets T n µ the role of sets T µ on which the density is non-zero.) We then use, for each given n the same scheme to cover the set X n with countably many T n µ , as was used in the proof of Theorem 1 to construct a countable covering of the set X ∞ , obtaining for each n a predictor ν n . Then the predictor ν is obtained as n∈N w n ν n , where the weights decrease subexponentially. The latter fact ensures that, although the weights depend on n, they still play no role asymptotically. The technically most involved part of the proof is to show that the sets T n µ in asymptotic have sufficiently large weights in those countable covers that we construct for each n. This is used to demonstrate the implication average KL divergence with probability 1. Of course, another very natural approach to the general problem posed above is to try and find predictors (in the form of algorithms) for some particular classes of processes which are of practical interest. Towards this end, the contribution of this work is in providing a specific form that some solution to this question has to have, if a solution exists: a Bayesian predictor whose prior is concentrated on a countable set. This is perhaps a rather simple form, which may be useful for constructing practical algorithms.

"if a set has a high µ probability then its ρ-probability does not decrease too fast, provided some regularity conditions." The proof is broken into the same steps as the (simpler) proof of Theorem 1, to make the analogy explicit and the proof more understandable.

Proof. Define the weights w k := wk -2 , where w is the normalizer 6/π 2 .

Step 1: densities. Define the sets

Using Markov's inequality, we derive

so that µ(T n µ ) → 1. (Note that if µ is singular with respect to ρ, as is typically the case, then ρ(x1..n) µ(x1..n) converges to 0 µ-a.e. and one can replace 1 n in (4) by 1, while still having µ(T n µ ) → 1.) Step 2n: a countable cover, time n. Fix an n ∈ N. Define m n 1 := max µ∈C ρ(T n µ ) (since X n are finite all suprema are reached). Find any

Observe that (for each n) there is only a finite number of positive m n k , since the set X n is finite; let K n be the largest index k such that m n k > 0. Let

As a result of this construction, for every n ∈ N every k ≤ K n and every x 1..n ∈ T n k using (4) we obtain

Step 2: the resulting predictor. Finally, define

where γ is the i.i.d. measure with equal probabilities of all x ∈ X (that is, γ(x 1..n ) = |X | -n for every n ∈ N and every x 1..n ∈ X n ). We will show that ν predicts every µ ∈ C, and then in the end of the proof (Step r) we will show how to replace γ by a combination of a countable set of elements of C (in fact, γ is just a regularizer which ensures that ν-probability of any word is never too close to 0).

Step 3: ν predicts every µ ∈ C. Fix any µ ∈ C. Introduce the parameters ε n µ ∈ (0, 1), n ∈ N, to be defined later, and let

, for any k > 1 and any n ∈ N, by definition of these sets. Since the sets

We can upper-bound µ(T n µ \T n j n µ ) as follows. First, observe that

Then, from (4) we get

Observe that for every n ∈ N and every set A ⊂ X n , using Jensen's inequality we can obtain

Thus, from ( 12) and ( 9) we get

where in the second inequality we have used the fact that entropy is maximized when all events are equiprobable, in the third one we used |X n \T n µ | ≤ |X | n , while the last inequality follows from (5). Combining (10) with the bounds (11), ( 13) and ( 14) we obtain

Since d n (µ, ρ) = o(n), we can define the parameters ε n µ in such a way that log ε n µ = o(n) while at the same time the bound (15) gives µ(T n µ \T n j n µ ) = o(1). Fix such a choice of ε n µ . Then, using µ(T n µ ) → 1, we can conclude

We proceed with the proof of

where the first inequality follows from (8), the second from ( 7), and in the equality we have used w j n µ = w/(j n µ ) 2 and j n µ = 1/ε µ n . Next we use the decomposition

From ( 17) we find

) where in the second inequality we have usedlog ε n µ = o(n) and d n (µ, ρ) = o(n), in the last inequality we have again used the fact that the entropy is maximized when all events are equiprobable, while the last equality follows from ( 16). Moreover, from (8) we find

where in the last inequality we have used γ(x 1..n ) = |X | -n and µ(x 1..n ) ≤ 1, and the last equality follows from (16). From ( 18), ( 19) and ( 20) we conclude 1 n d n (ν, µ) → 0.

Step r: the regularizer γ. It remains to show that the i.i.d. regularlizer γ in the definition of ν (8), can be replaced by a convex combination of a countably many elements from C. Indeed, for each n ∈ N, denote

for each x ′ 1..n ∈ A n , n ∈ N, and let γ ′ := k∈N w k γ ′ k . For every µ ∈ C we have

for every n ∈ N and every x 1..n ∈ A n , which clearly suffices to establish the bound II = o(n) as in (20).

Discussion

For two measures of quality of prediction that we have considered, namely, total variation distance and expected average KL divergence, we have shown that if a prediction for a class C of measures exists, then a predictor can also be obtained as a Bayesian mixture over a countable subset of C. The first possible extension of these results that comes to mind is to find out whether the same holds for other measures of performance, such as prediction in KL divergence without time-averaging, or with probability 1 rather then in expectation. Maybe the same results can be obtained in more general formulations, such as f -divergences of [Csiszar, 1967]. More generally, the questions we addressed in this work are a part of a larger problem: given an arbitrary class C of stochastic processes, find the best predictor for it. One can approach this problem from other sides. For example, the first question one may wish to address is for which classes of processes a predictor exists; see [Ryabko, 2008] for some sufficient conditions, such as separability of the class C. Another approach is to identify the conditions which two measures µ and ρ have to satisfy in order for ρ to predict µ. For prediction in total variation such conditions have been identified [Blackwell andDubins, 1962, Kalai andLehrer, 1994] and, in particular, in the context of the present work, they turn out to be very useful. [START_REF] Kalai | Weak and strong merging of opinions[END_REF]] also provides some characterization for the case of a weaker notion of prediction: difference between conditional probabilities of the next (several) outcomes (weak merging of opinions). In [START_REF] Ryabko | Predicting non-stationary processes[END_REF]] some sufficient conditions are found for the case of prediction in expected average KL divergence, and prediction in


Keywords: R&D investment, Technology Spillovers, Complementarities, Economic growth, Dynamic Panel Data, Cointegration, Unit root test, Private investment, Public investment, R&D cooperation

This paper investigates the relationship between private and public investment in R&D, while taking into account the effect of several instruments policies such as subsidies and taxes. We design a new look of knowledge spillovers and R&D cooperation to explain the contribution of public and private R&D on growth. We propose a heterogeneous dynamic panel data model to consider the endogenous effect of R&D investment. We also distinguish between the estimated long and short run results. Our results based on a sample of 23 countries over the period 1992-2004 indicate that both public and private investments in R&D are complementary. By establishing an endogenous growth model, the estimates indicate that public and private R&D depends on the host country's human capital investment. Results indicate that foreign direct investment is a more significant spillover channel than imports.

I. INTRODUCTION

Is public R&D complementary to private R&D, or does it substitute for and tend to "crowd out" private R&D? Conflicting answers are given to this question. We survey the body of available econometric evidence accumulated over the past 16 years. A framework for analysis of the problem is developed to help organize and summarize the findings of econometric studies based on dynamic panel data from various countries 1 . We conclude by offering suggestions for improving future empirical research on this issue.

Most people think that government R&D activities contribute to innovation and productivity, many economists and policymakers have grown frustrated with the paucity of systematic statistical evidence documenting a direct contribution from public R&D (see Paul et al,. [START_REF] Paul | Is public R&D a complement or substitute for private R&D? A review of the econometric evidence[END_REF]. Econometric findings concerning the productivity growth effects of R&D seems to be that there is a significantly positive and relatively high rate of return to R&D investments at both the private and social levels (Cassiman and Veugelers, [START_REF] Cassiman | In search of complementarity in Innovation Strategy: Internal R&D, Cooperation in R&D and external technology Acquisition[END_REF]). In a recent survey, Paul et al,. [START_REF] Paul | Is public R&D a complement or substitute for private R&D? A review of the econometric evidence[END_REF] suggest that the especially pronounced differential over the returns on tangible capital investments observed at the private level may reflect individual firms' perceptions of especially high private risk in the case of R&D. Public funding of R&D can contribute indirectly, by complementing and hence stimulating private R&D expenditures. Our approach will be to adopt a new econometric approach using a dynamic panel data studies to analyse if public investment in R&D are complement or substitute for private investment in R&D. In literature review, we can conclude that the majority of the econometric studies are concentrated on the impact of public R&D contracts and grants upon private R&D investment by manufacturing firms and industries (see for example, Lach [START_REF] Lach | Do R&D subsidies stimulate or displace private R&D? Evidence from Israel[END_REF], Christopher [START_REF] Christopher | Market structure and endogenous productivity growth: how do R&D subsidies affect market structure?[END_REF] and Eric [START_REF] Eric | R&D efficiency and economic performance: A cross-country analysis using the stochastic frontier approach[END_REF]).

The object of our paper is to give the theoretical and empirical arguments which allow a satisfactory apprehension of the role that the authorities must play in the fields of research and innovation. The activity of R&D represents a significant source of development of new knowledge and technological innovation (Guellec and Van-Pottelsborghe, [START_REF] Guellec | Does Government Support Stimulate Private R&D?[END_REF]). The effort towards activity of R&D involves with a great importance and this through several resources devoted to the various sectors and institutions of research. Expenditure of research and development especially constitute a principal source of growth of productivity for innovating countries. Whereas, Sigrid [START_REF] Lopez-Bazo | Complementarity between local knowledge and internationalisation in regional technological progress[END_REF] and Ting [START_REF] Ting | Foreign direct investment and growth under economic integration[END_REF] suggest that, for countries, where the activity of R&D misses almost technological knowledge and innovations of which they profit are generally resulting from the importation of equipment and goods of intensives investments in technical progress. At the same times, Chaturvedi and Chataway [START_REF] Chaturvedi | Strategic integration of knowledge in Indian pharmaceutical firms: creating competencies for innovation[END_REF] recommend that knowledge capabilities and knowledge management can be considered as key resources for firms in both developed and developing countries.

II. R&D INVESTMENT AND THE GOVERNMENT'S R&D POLICIES

Today, we can observe an expansion of policies of innovations in the developed countries which devote great investment for R&D. What proves the creation of the climates favourable to the level of these countries for the innovation? It is significant that during these last years, companies of high technology or advanced technology's (pharmaceutical, aeronautical…) expenditure of research and development increased significantly. The role of the governmental policies as regards R&D is not to neglect. Indeed, the policies of innovation define specific actions of the State, which must encourage the accumulation of a qualified labour on the one hand, and to help the companies to prospect the markets on the other hand. This justifies the need for the public administrations for supporting the R&D.

Thus, which are the reasons of the government aid and the mechanisms the alternate ones available to the public administrations to support the R&D? To answer these questions we try to analyze the justification of the government aid with the R&D starting from the economic theories of growth.

Neoclassic growth theory

For Neo-classic theory of growth, technical progress is supposed to be exogenous factors. With the balance of long term, population growth and technical progress determine the level of the growth rate. This implies, according to the basic assumptions, that the long-term growth rate is stable, and given in an exogenous way. Within this framework, the impact of an action of the authorities is practically ignored.

The neoclassic theory of the growth supposes that the economy starts from a weak relationship between capital and labour. Just as the marginal returns on capital are decreasing. What reduces the Complementarity between private and public investment in R&D: A Dynamic Panel Data analysis Tarek Sadraoui* & Naceur Ben Zina encouragement to be invested in the new capital? Thus each new unit of capital produces a lower income and less large savings. In the long run, there will be absence of incentive to invest. In short, we can say that the assumptions which underlie the neo-classic theory are not realistic. The technological change is not always an exogenous factor outside the market, determined by an unknown process. To the 20th centuries, a good number of discoveries and progress were carried out in the commercial sector by companies with lucrative goal and not by public administrations or universities where research is directed by non-commercial forces. Markets are seldom in perfect competition, moreover, the private sector is not capable to produce all the desired goods and services, because some of them are goods public and certain others produce external effects.

Endogenous growth theory

The endogenous theory of growth recommends the relaxation of certain neo-classic assumptions and incorporates the failures of the market. However, the economic growth in the long run is directed by the accumulation of the factors of production founded this faith on knowledge, in particular, human capital, training, R&D and innovation (Griliches, [START_REF] Griliches | R&D and Productivity: Econometric Results and Measurement Issues[END_REF]). The endogenous models of growth are characterized by a great diversity of the resources selected: The investment in physical capital, in human capital, public capital, and labour division, learning by doing, research and the technological innovation (Romer,[START_REF] Romer | Endogenous technological change[END_REF]).

R&D investment and market imperfection

Economic theory and empirical proof show that technical progress, because of its incidence on the factors of production, constitutes key element in the long run determining economic growth; in certain countries, it represents even the most significant element. However, it is not a question of an economic justification of the official intervention for allocate the resources in favour of R&D. But, this intervention in a market economy is justified by incapacity of market to distribute resources in an efficient or acceptable way as regards social aspects. With regard to the investment in R&D, external effects and market imperfections testify the incapacity of market, and the effects are felt not only beyond particular companies but also beyond national borders.

In a market economy, a company will not invest in a project if it knows that it can not adapt the possible receipts, however if it cannot adapt a portion of these receipts, it will invest if this portion is enough to make a profitable investment. Asymmetrical information and imperfect competition constitute two other kinds of imperfections of market involving under investment in R&D. For example, asymmetrical information prevents effective operation of capital market. Indeed, it can involve rationing of appropriations as well as abandonment of investments in R&D in projects with strong chances of success thanks to the plan of financing, and the continuation of investments in the project having weak chances of success.

III. RESEARCH METHODOLOGY

Our study contributes to the empirical literature -which is discussed here-on the analysis of the existence of a relation between private and public investment in R&D and their real effect on economic growth; do public funds substitute or complement private R&D expenditure?

We derive our econometric specification from a function including interactions between internal and external R&D in the augmentation of the knowledge stock. The model also takes into account potential productivity convergence by including lagged productivity levels. Our study's inferences are based on a dynamic panel data model, which allows us to control for the existence of unobserved fixed effects that are likely to affect R&D decisions. Estimation is carried out by several consistent dynamic panel data methods, among which generalized method of moments, which allows for the presence of weakly endogenous explanatory variables. In this way the analysis can take into account both degree and possibility effects of R&D to address the issue of optimal combinations of R&D expenditures.

In this paper, we contribute the first panel data study exploring complementarity between public R&D and private R&D in a dynamic panel framework. We examine the impact of internal and external R&D on economic growth in sixteen-year panel for 23 developed and developing countries.

IV. COMPLEMENTARITY VERSUS SUBSTITUTABILITY BETWEEN PRIVATE AND PUBLIC R&D

Theoretical work did not succeed in slicing on favourable or unfavourable effect using certain political instruments on the level of R&D in private sectors. The results of each model strongly depend on its structure and its assumptions. Empirical work, leads to homogeneous results and identifies a positive effect of public R&D on that private (Paul, et al,. [START_REF] Paul | Is public R&D a complement or substitute for private R&D? A review of the econometric evidence[END_REF]). With an aim of knowing the relation between public and private R&D we give an overall picture of the activities of R&D in world. Indeed, in this section, we attach more importance to activity of public and private R&D in the most significant poles in world. Through time and with the improved scientific methods in particularly studies of Jason, et al., [START_REF] Jason | A model of R&D valuation and the design of research incentives insurance[END_REF], it became clearly that the final situation towards the effect of the public funds of R&D cannot be made. Thus, in general, two fields can be identified and which are used to analyze the relation between private investment and public in research and development with knowing quantitative and qualitative studies: On the one hand; for the qualitative studies, data are frequently based on the investigations. On the other hand; for the quantitative studies, they are based on macro and micro-economic information of a significant number of companies (Cassiman and Veugelers, [START_REF] Cassiman | In search of complementarity in Innovation Strategy: Internal R&D, Cooperation in R&D and external technology Acquisition[END_REF]).

Today, several activities of R&D are carried out on the level of the services sector. On the one hand, this is due to the external sources of the strategies of manufacturing industries in the Eighties. On the other hand, the transformation of information and technology of communication get more opportunities for innovating sectors. So the governments help more and more activities of R&D in several sectors with an aim of stimulating technological performances of their countries. Thus, several examples can be quoted. At this level, for the Nineties and more precisely in 1999, the total expenditure of R&D of Germany is 47 billion dollars where 66% of this amount is invested by private industries, 18% by government and the remainder are invested by foreign companies.

Thus, Claudia, [START_REF] Claudia | Complementarities of innovation activities: An empirical analysis of the German manufacturing sector[END_REF] suggest that an international comparison on behalf of public programs of R&D shows that Germany is one of principal countries which grant funds for the technological performance. At this level, manufacturing industry plays a very significant role concerning R&D. For example, the strategic planning of the national research evaluation in Thailand as indicates it Jarunee, [START_REF] Jarunee | Building effective research evaluation systems to assist R&D investment decisions[END_REF] is to allocate the budget to support the research programmes and projects. Jarunee suggest that to improve the model evaluation framework for R&D investments, the public hearing forum was organised. From there, a question emerges up to what point evolution of public funds of R&D makes it possible to stimulate R&D carried out by private sector, and on which level results are checked? Recently, an econometric micro study tackled the question of the impact of political instruments about activity of R&D deprived on the level of companies.

Several other studies are more precisely interested in testing the effects of public subsidies in R&D on the amount of deprived investment see for example Lai et al., [START_REF] Lai | Technology spillovers, absorptive capacity and economic growth[END_REF], Hans and Almas, [START_REF] Hans | The impact of public funds on private R&D investment: New evidence from a firm level innovation study[END_REF] and Christopher, [START_REF] Christopher | Market structure and endogenous productivity growth: how do R&D subsidies affect market structure?[END_REF]. The major goal of these studies is to know if public subsidies of research and development can have an effect of reduction or increase in the expenditure of R&D. Most of results suggest that public subsidies of R&D on the level of several industries showed that there is a small tendency to the effect of ousting "Crowding out". In addition, it seems not to have any effect or degree of Complementarity.

In the following section we empirically test fundamental relation which we seek to analyze in the case of 23 countries for the period 1992-2004, in other words we test the existence of a relation of Complementarity and to check this result.

V. EMPIRICAL VALIDATION: DYNAMIC PANEL DATA

The objective of this work is to test the impact of an action of public policies empirically on the evolution of R&D in private sector while trying to surmount limits. The modelling that we follow to measure the effect of the R&D deprived on the public one; while taking into account some determinants of private R&D; is the one of Bettina and al,. (2002). This modelling has also been applied by: Busom, [START_REF] Busom | An empirical evaluation of the effects of R&D subsidies[END_REF] and Lach, [START_REF] Lach | Do R&D subsidies stimulate or displace private R&D? Evidence from Israel[END_REF]. The gait of these authors can be summarized like follows:

Private R&D = ß * public R&D + control variables + e

The underlying logic is simple: If the coefficient β* has a positive sign we can say that public R&D are complementary for private R&D. In other words, an increase of 1% of public research and development level entails a growth of β*% of private R&D. On contrary, if β* has a negative sign we can say that there is a relation of substitutability between public and private R&D. In this part we try, to give a general setting for the models to estimate while putting accent on some remarks and inconveniences of these models. We apply a dynamic panel data model. Finally, after having estimated the model we analyze results.

Dynamic panel data: Definition and evaluations method

Dynamic models are characterized by presence of one or several endogenous variables delayed among explanatory variables. Our specified model is a dynamic panel model is given by:

' -1    y y x v it it it it   (1)
Under another forms one was writing our model as below:

-

1 1 2 3 4        G V A R R M IDE it it it it it it it i        (2)
Where;

-1 y it : Endogenous variable appears in the regression as being a retarded explanatory variable. In other words, present stocks of research and development of country (i) are explained by stocks of research of the period (t-1). X: Represent the vector of exogenous variables; these variables are added value (VA), public research (G), import (M), foreign direct investment (IDE) and private research;

(α, β): Designate parameters to estimate; μ i : Constitute individual heterogeneity as: μ i i.i.d. ~ N [(0, 1)]; And: v i,t is stochastic term as: v i,t ~ i.i.d. [(0, 1)]. y i,t is the logarithm of volume of R&D in country (i).

x i,t is determinant vector of R&D. μ i is the specific effect of country (i). This specific effect can be a stationary or uncertain effect.

Evaluation Method

The evaluation of the model by traditional methods (Ordinary Least Square "OLS" and within) gives biased and non convergent values because of inter-relationship between retarded endogenous variable and individual heterogeneity. We try to demonstrate for the case of a simple model the inconveniences of these methods of evaluations.

For dynamic panel model, Within transformations and Ordinary Least Squares are biased and non-convergent estimators.

We assume the simple specification:

-1 y y it i it it       Let us pose 1
  and that y it is stationary then: 1) )

1 1 1 1 1 1 ˆ2 T N i it it t T N i it t y y y a = - = = - = = = å å å å               N i T t it N i T t it y y v it i 1 1 1 1 1 1 2    ˆ is convergent if   0 1 lim 1 1 1           N i T t it N y v it i NT p             N i T t it N y v it i NT p 1 1 1 1 lim  =             T i i T T T T T y        ) 1 ( ) 1 ( , cov 1 1 1 2 0 2 1 2 l i m 1 1 1 T N p y i i t N T t N           2 2 2 2 1 1 1 1 1 0 1 2 2 2 1 1 (1 ) 1 N t t t y i i T T N T T                              2 2 2 1 1 c o v ( , ) 0 2 (1 ) 1 1 t t y i i T                       2 2 2 ( 1) 2 2 ((
T T T y T               
In summary, the bias is positive and increases with the variance of the specific effect. Indeed, y i,t is function of v i,t and y i,t-1 is also. y i,t-1 is an explanatory variable correlated with stochastic term. It introduces a bias in the value of ordinary least squares. Even as putting hypothesis that stochastic terms are not correlated, this value is non-convergent.

For within case we consider the following transformation:

    - - -. -1 . -1 .-1 y y y y it i it i it i      With:     0 s i i j e t si t s 0 E v v it js E ifn o t v v it js           
While posing as in Baltagi [START_REF] Baltagi | A Dynamic Demand Model for Liquor: The Case for Pooling[END_REF]:

1 . 1 1
y i t y i T     
Thus we can write that:

     . 1 . 1 1 1 ˆ2 1 . 1 1 1 T N y y y y it i it i i t T N y y it i t i                        . 1 . 1 1 1 2 1 . 1 1 1 T y y it i it i N i NT t T y y it i N i NT t                    
The numerator is convergent when the second term converges towards zero. The numerator of the second term:

               N i T t i it N y y i it NT p 1 1 . 1 . 1 1 lim   2 2 2 ( 1)
( 1 )

T v T T T a s a a - - - + = - And the denominator             N i T t i it N y y NT p 1 1 2 1 . 1 1 lim 2 2 2 2 ( 1) 1 2 1 1 (1 ) T v T T T T a a s a a a ae ö ÷ - ç - - + ÷ ç ÷ ç = - ÷ ç ÷ ç - ÷ ÷ ç ç ÷ è ø -
The estimator of Least Square Dummy variable is convergent if T is infinite. If T is fixed, N, the estimator is non-convergent. Our model should not be estimated by the method of OLS and LSDV due to the fact that estimating by these methods led to ad hoc results. Which are then adequate methods to estimate our model? We propose below two methods which consist in obtaining consistent estimators.

Anderson and Hsiao Method

Anderson and Hsiao [START_REF] Anderson | Formulation and Estimation of Dynamic Models Using Panel Data[END_REF] proposed, initially, to write the model from first difference to eliminate individual heterogeneity. They propose for the transformation two instruments.

      - - -1 -2 .-3 1 3 (1) (3) ˆ- - .-2 -2 .-3 1 1          T N y y y y it it it i i t vi T N y y y y it i it i i t  And     - -1 -2 1 3 ˆ( 2 ) -.-2 -2 1 1          T N y y y it it it i t v i T N y y y it i it i t  (4)
The two values are convergent when N and T   .

However, an inter-relationship always persists between endogenous variable in first difference and residual term. Authors proposed to resort to the method of instrumental variables to surmount this problem. Thus, they propose to use instrument endogenous variable with two lags or his first differences. These instruments are correlated with explanatory variable and are not with residual term. To get more efficient results, Arellano and Bond [START_REF] Arellano | Some test of specification for panel data: Monte Carlo evidence and application to employment equations[END_REF] approach permits to get a value of generalized moments "GMM" more efficient.

Arellano and Bond Approach

Arellano and Bond [START_REF] Arellano | Some test of specification for panel data: Monte Carlo evidence and application to employment equations[END_REF] are the first in 1991 that proposed an extension of GMM introduced initially by Hansen [START_REF] Hansen | Large Sample Properties of Generalized Method of Moments Estimators[END_REF], to the case of panel data for a simple model AR (1):

(5) -1    y y v it i it it   Where 0  
We consider the case where temporal dimension is small while individual dimension (N) is important. We consider that individual effects are stationary and we assume traditional hypotheses of residues:

In difference models (5) can be written as below:

0 (6) -1      y y u it it it   Where u it = v it -v it-1 .
We test for every individual of the linear restrictions of type:

  - 0 2,..., ; 3,... -1 -            E f o rj t t T y y y it it it j  (7)
The gait of Arellano and Bond, in presence of the exogenous variables, consists in estimating the model in difference:

'( ) (8)
( -) 1         p y v X L y it it it k i t k k  
Moment conditions and instruments matrix are given respectively by:

    0 2, 2, 3, .... 0 2, 2, 3,....                  is it it it E p o u r t T E p o u r t T y v v X     (9) 1 2 3 0 0 . . . . 0 0 . . . . . . . . . . . . . . .. . . . .                      ip ip iT ip ip Z Z Z X X X (10)
The preceding dynamic model ( 8) can be rewritten for each individual in the following form:

(11)    y W V i i i i i   
Where τ is a vector of parameter and W i is a matrix that contains the retarded dependent variable and explanatory variables. The method proposed by these author's permits to get a GMM in two stages is written in following form:

-1 * * *' *' ˆ(12) '                                             W y Z Z W W Z A Z A i i i N i N i i i i i i i i
 However, to have previous value GMM, it is necessary to pass by a first stage that consists in making wished transformation (first difference or orthogonal deviation), to find and to use instruments matrix and to achieve a first evaluation named "evaluation of first stage". This stage corresponds to an evaluation that permits to provide estimated residues after transformation. In the first stage, the values are gotten while using H i as: * *' ˆˆ(13)  v v H i i i

               H i And   -1 1 ' ' ( 1 4
)

           A Z Ä H Z N Z H Z I i i i N N i
The objective of transformation is, as at Anderson and Hsiao [START_REF] Anderson | Formulation and Estimation of Dynamic Models Using Panel Data[END_REF], to eliminate individual heterogeneity of the model. The number of instrument increases in the time for every individual. In the case where exist explanatory variables x it in the model correlated with heterogeneity individual µ i . Optimal instruments matrix corresponding Z i is equal to:

               y x x i i i y x x x i i i i y y x x i T i i i T -1) 1 1 ( -2)
Arellano and Bond [START_REF] Arellano | Some test of specification for panel data: Monte Carlo evidence and application to employment equations[END_REF] propose a test verifying the absence of autocorrelation of first and second order. Thus, if distribution is non auto-correlated, this test gives a value of residues differentiated negative and significant to first order and non significant to the second order. This test that is based on auto-covariance of residues follows a normal law N (0,1) under hypothesis H 0 . Otherwise, authors propose the test of validity of instruments of Sargan (1988a). The statistical test equal to:

* * ' ˆ(1 5 ) ˆ'                    A v S N Z v Z i i i i

VI. EMPIRICAL RESULTS

The unit root tests became a current step for analysis of time series stationnarity. However, practical application of these tests on panel data is recent. The tests most frequently used are those of Levin and Lin, [START_REF] Levin | Unit root tests in panel data: asymptotic and finite sample properties[END_REF] and of Im, Pesaran and Shin [START_REF] Im | Testing for unit roots in heterogeneous panels[END_REF] 3 . Recently, several procedures of unit root tests and Cointegration were developed for panel data models. The addition of individual dimension to temporal dimension offers an advantage, in practical application of unit root and Cointegration tests (Pedroni,[START_REF] Pedroni | Critical Values for Cointegration Tests in Heterogeneous Panels with Multiple Regressors[END_REF][START_REF] Pedroni | Panel Cointegration: Asymptotic and Finite Sample Properties of Pooled Time Series Tests with an Application to the PPP Hypothesis[END_REF]).

The checking of non-stationary properties for all panel variables leads us to study the existence of a long run relation between these variables. The Cointegration study by applying Pedroni Cointegration tests based on unit root tests on residues estimated. Cointegration tests on panel data consist in testing the presence of unit root in the estimated residues. However, the problem of fallacious regressions, of the time series, also arises in the case of panel data. First step is to test unit root for each of series.

Unit Root Tests

Levin and Lin [START_REF] Levin | Unit root tests in panel data: asymptotic and finite sample properties[END_REF], consider the following model:

y i,t =  i y i,t-1 + Zit  + ui,t (i=1, …, N; t=1, …, T) (16)
Where, Z i,t is the deterministic component and u i,t is a stationary process.  i is the fixed effect, The Levin and Lin, [START_REF] Levin | Unit root tests in panel data: asymptotic and finite sample properties[END_REF] tests assume that ui,t are iid (0, 2 u ) and i= for all i. The LL test is restrictive in the sense that it requires  to be homogeneous across i. Im, Pesaran and Shin [START_REF] Im | Testing for unit roots in heterogeneous panels[END_REF] (IPS) allow for a heterogeneous coefficient of y i,t-1 and propose an alternative testing procedure based on averaging individual unit root test statistics. IPS suggested an average of the augmented Dickey-Fuller (ADF) tests when u i,t is serially correlated with different series. Correlation properties across cross-sectional units, i.e;

, 1      i p i t j ij it j it u u   .
Substituting this u i,t in (1) we get:

, 1 1          p i y y y z i t i it j ij it j it it     (17)
The null and for all countries i the alternative hypothesis are: Ho:

 i = 1 Ha:  i < 1
For at least one i, the IPS t-bar statistic is defined as the average of the individual ADF statistic as:

i N i t N t  1 1   
Where ti is the individual t-statistic of testing Ho: i = 1 in [START_REF] Jarunee | Building effective research evaluation systems to assist R&D investment decisions[END_REF]. It is known for a fixed N as T  

1 0 1/ 2 1 2 0 iz iz i i T iz W dW t t W            (18)
IPS assumes that t iT are iid 4 are have finite mean variance. Then;

1 [ / 1 1 (0, 1) [ / ] 1               N N t Et iT iT i N i N Var t iT i   (19) As N  central limit theorem. Hence   [ / 1 (0, 1) [ / ] 1      N t E t iT i t N IPS Var t iT i   (20)
As T  followed by N  sequentially, the values of E[tiT/i=1] and Var[tiT/i=1] have been computed by IPS simulations for different values of T and is. As applying test on our complete model our results is summarized in table 1

VII. CONCLUSION

In our survey, we tried to put accent on private and public investment in R&D, for the case of 23 countries which presents different levels of R&D. We tried to clarify relation that exists between private and public research. This empirical survey wanted to give account, the effects of different determinants on private investment in R&D and to know if public and private investments in R&D are complement or substitute.

Econometric approach consists in the regression of some measures of private R&D on public R&D with some control variables. The evaluation that we presented in our work corresponds to GMM evaluation in first difference and in orthogonal deviation. We prefer to refer to results of this evaluation because it permits to eliminate rigorous way all bias to none observed individual heterogeneity and offer, a better efficiency of results. Empiric evaluations confirm a positive effect of public R&D in different country (positive and meaningful effect in all evaluations). Results of our empiric survey are relative for our sample and they go in the sense of results of ulterior studies, which showed that there is a positive and meaningful relation between private and public investment in R&D.

All results are in favor of a positive relation between private and public R&D which can be assumed by a complementarities between them. In our study we have indicate that all variables are stationary by an application of unit root tests that's can contribute to search Cointegration relation between them and determinate the number of these relation. Another important think, we can give the impact of public R&D to private R&D for each country to specify the nature of relation and how private R&D contributes for public sector. In summary, all countries must stimulate private sector in R&D activities to promote economic growth and integrate a new innovation system which can go with their own economic environment.

Some studies put in value of other factors that can be important as: competition in the market, public politics and cooperation concerning R&D between firms. Cooperation in R&D is a part of the new strategies developed by firms in more global and competitive economic environment. These last factors are not to disregard and can be subject of a future research concerning the relation between public and private investment in R&D.

APPENDIX

  ...........-1 2

Table 1Unit

 1Unit

		root Tests
	Statistics	R	G	M	VA	IDE
	Levin-Lin ADF stat	-2.357	-0.343	1. 809	1.489	-1.674
	IPS ADF stat	-2.296	1.773	2.175	1.640	-1.572


Claire Bost

Georges-Henri Cottet

Emmanuel Maitre

Linear stability analysis of a Level Set model of immersed elastic membrane

. We also present numerical tests on the full model which further illustrate the relevance of this simplified model.

Introduction

We study the numerical stability of an elastic membrane immersed in an incompressible 2D fluid. This coupled model is described by an Eulerian version of the Immersed Boundary Method of Peskin [START_REF] Peskin | The immersed boundary method[END_REF], proposed by Cottet & Maitre ([7], [START_REF] Cottet | A Level Set method for fluid-structure interactions with immersed surfaces[END_REF]), and using a Level Set formulation. Indeed, the stability of the Immersed Boundary Method is well known to be affected by the application of the elastic force on the structure, as soon as it is treated explicitely. In order to remove the time step restrictions, many works have been dedicated to the study of implicit schemes ( [START_REF] Tu | Stability and instability in the computation of flows with moving immersed boundaries: a comparison of three methods[END_REF], [START_REF] Mayo | An implicit numerical method for fluid dynamics problems with immersed elastic boundaries[END_REF], [START_REF] Stockie | Stability Analysis for the Immersed Fiber Problem[END_REF], [START_REF] Stockie | Analysis of Stiffness in the Immersed Boundary Method and Implications for Time-stepping Schemes[END_REF], [START_REF] Newren | Enhancing the Immersed Boundary Method: stability, volume conservation, and implicit solvers[END_REF]). However these schemes seem to be unusable in practice because of the large computational cost related to the iterative resolution of the strongly non-linear coupling at each time step. The definition of semi-implicit or approximative implicit schemes allows to obtain more realistic methods. See [START_REF] Tu | Stability and instability in the computation of flows with moving immersed boundaries: a comparison of three methods[END_REF], [START_REF] Stockie | Stability Analysis for the Immersed Fiber Problem[END_REF], [START_REF] Stockie | Analysis of Stiffness in the Immersed Boundary Method and Implications for Time-stepping Schemes[END_REF] [START_REF] Fernandez | A projection semiimplicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF] et [START_REF] Boffi | Numerical stability of the finite element immersed boundary method[END_REF] and the references therein for details about such schemes. There exists several ways to carry out the study of numerical stability. Among them is the control of discrete energy conservation ( [START_REF] Boffi | Numerical stability of the finite element immersed boundary method[END_REF], [START_REF] Newren | Enhancing the Immersed Boundary Method: stability, volume conservation, and implicit solvers[END_REF]), which is based on the corresponding physical property. Another way is to perform a linear analysis near an equilibrium state of the system, in order to identify the modes associated to the membrane. In ( [START_REF] Stockie | Stability Analysis for the Immersed Fiber Problem[END_REF], [START_REF] Stockie | Analysis of Stiffness in the Immersed Boundary Method and Implications for Time-stepping Schemes[END_REF]) this approach enables the authors to show that the instability is in particular increased by the combination of a small fluid viscosity and a strong elastic force. The aim of this work is the derivation of stability conditions on the time step for schemes which treat the coupling explicitely, semi-implicitely or implicitely. We perform the analysis on a 1D linearized version of the model as its original form is not very suitable for numerical analysis. More precisely, the elastic force on the membrane is reduced to a surface tension. Despite of these simplifications, this ad-hoc model enables us to exhibit some clear relations between the fluid viscosity, the elastic force, and the numerical parameters. The results are in good accordance with some conditions of the litterature, in particular those proposed either by Brackbill et al. [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF] or Vigneaux [START_REF] Vigneaux | Méthodes Level Set pour des problèmes d'interface en microfluidique[END_REF], depending on the physical parameters range. In section 2 we briefly introduce the 2D Eulerian Immersed Boundary model and we derive the 1D linearized one used for the analysis. The time-stepping schemes studied are also presented. Section 3 is then dedicated to the computation of necessary stability conditions and to the comparison of these conditions with existing ones. Finally, numerical experiments on a relaxation test validate in section 4 the analytical results in the non-linear 2D case, for some range of parameters.

Immersed elastic membrane problem

Let Ω be an open bounded domain of R 2 , filled with a viscous incompressible and homogeneous fluid of density ρ > 0 and viscosity µ > 0. Inside this domain, we consider an immersed elastic membrane Γ e (t) during a time interval [0, T ], T > 0. The membrane is supposed massless, and external forces are neglected. Let u and p be the fluid velocity and pressure fields. We introduce the Level Set function φ initialized as the signed distance to Γ e (0), and such that :

Γ e (t) = {x ∈ Ω, φ(t, x) = 0}, ∀t ∈ [0, T ].

Eulerian Immersed Boundary model in dimension 2

We study the following Eulerian version of the Immersed Boundary Method for this model [START_REF] Cottet | A Level Set formulation of immersed boundary methods for fluid-structure interaction problems[END_REF], [START_REF] Cottet | A Level Set method for fluid-structure interactions with immersed surfaces[END_REF]:

             ρ (u t + u.∇u) + ∇p -µ∆u = F e F e = ∇(E ′ e (|∇φ|)) • ∇ × φ |∇φ| • ∇ × φ |∇φ| -E ′ e (|∇φ|)κ(φ) ∇φ |∇φ| |∇φ| ε ζ φ ε div u = 0 φ t + u • ∇φ = 0
(1) with the linear stress-strain relationship

E ′ e (r) = ν e (r -1) (2)
where ν e and κ(φ) are the stiffness coefficient and the mean curvature of Γ e (t), ε is a numerical parameter related to interface smoothing, ζ is a cut-off function,

and F e is finally the elastic force arising from the presence of the membrane in the fluid. For a complete description this system has to be supplemented with initial and boundary conditions. One common choice is periodic or homogeneous boundary condition for the velocity, so that no boundary condition is needed for φ.

Remark 2.1. When using this formulation, particular care must be taken of the spatial discretization of the elastic force term. Indeed, some derivatives have to be developed in order to avoid strong numerical instabilities, which must be kept apart from the instabilities related to the time-stepping choice. The study of these numerical instabilities is beyond the scope of this work and has been developed in [START_REF] Bost | Méthode Level Set et pénalisation pour le calcul d'interactions fluide-structure[END_REF].

1D linearized model

We now introduce a 1D linearized version of the previous model. The elastic force is first reduced to a surface tension by taking E ′ (r) = ν e in (1). This results in the following 2D system:

         ∂u ∂t + u.∇u -µ∆u + ∇p = -ν e κ(φ)∇φ 1 ε ζ φ ε , div u = 0, ∂φ ∂t + u.∇φ = 0. (3)
which next gives rise to the 1D model:

     ∂u ∂t + u ∂u ∂x -µ ∂ 2 u ∂x 2 = -ν e ∂ 2 φ ∂x 2 ∂φ ∂x 1 ε ζ φ ε , ∂φ ∂t + u ∂φ ∂x = 0. (4)
We investigate the stability of (4) around the trivial stationary solution (u, φ) given by: u

(x) = 0, φ(x) = x. (5)
Let (u, φ) = (u + u, φ + φ) be a solution of (4), where ( u, φ) is a small perturbation. Plugging these expressions in (4) we get:

         ∂ u ∂t + u ∂ u ∂x -µ ∂ 2 u ∂x 2 = -ν e ∂ 2 φ ∂x 2 1 + ∂ φ ∂x 1 ε ζ φ + φ ε , ∂ φ ∂t + u ∂ φ ∂x + ũ = 0. (6)
By getting rid of the non-linear terms, we then obtain:

       ∂ u ∂t -µ ∂ 2 u ∂x 2 = -ν e ∂ 2 φ ∂x 2 1 ε ζ φ ε + o( u, φ), ∂ φ ∂t + u = o( u, φ). (7)
Actually,

∂ 2 φ ∂x 2 1 + ∂ φ ∂x 1 ε ζ φ + φ ε = ∂ 2 φ ∂x 2 1 ε ζ φ + φ ε + o( φ) = ∂ 2 φ ∂x 2 1 ε ζ φ ε + φ ε ζ ′ φ ε + o( φ) = ∂ 2 φ ∂x 2 1 ε ζ φ ε + o( φ).
The remaining simplification concerns the function ζ. We define it by:

   ζ(r) = 1 if |r| ≤ 1 2 , ζ(r) = 0 else (8)
Let us point out that for ζ φ ε = 0, namely far from the fluid-membrane interface, the system (7) is no longer coupled. For the study below, we thus suppose that ζ φ ε = 1 everywhere, which means that we limit the study in a neighbourhood of the membrane. We are clearly conscious of the restriction brought by this hypothesis, since it is not clear how the presence of the membrane affects the stability in the whole computational domain, and conversely whether the computations in the remainder of the domain might bring stabilizing effects.

Finally, the 1D linearized model used for the stability analysis reads (we remove from now on the tilde on unknowns):

         ∂u ∂t -µ ∂ 2 u ∂x 2 = - ν e ε ∂ 2 φ ∂x 2 on [0, T ] × R ∂φ ∂t + u = 0 on [0, T ] × R u(0, x) = f (x), φ(0, x) = g(x) on R (9)
where ε > 0, µ ≥ 0 and ν e > 0 (we rule out the case ν e = 0 in order to keep [START_REF] Fernandez | A projection semiimplicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF] in a coupled form). This quite simplified model is not supposed to describe in a general way the fluid-membrane interactions, as the non-linearities, including in particular the inertial effects, have been neglected. However it enables us to perform a stability analysis and to draw some attractive conclusions.

Time-stepping schemes

Let ∆t > 0 and ∆x > 0 be the time and space discretization steps, x j = j∆x, j ∈ Z, and t n = n∆t, n ∈ N. We denote by u n j ≈ u(t n , x j ) and φ n j ≈ φ(t n , x j ), j = 0, ..., M + 1, the numerical approximation of solutions of (9), by some finite difference scheme to be described. Discretizing diffusion implicitely, so that the stability analysis will focus on the coupling, we define the following schemes to compute these numerical values:

• Scheme 1: Explicit (u, φ)-coupling:            u n+1 j -u n j ∆t -µ u n+1 j+1 -2u n+1 j + u n+1 j-1 (∆x) 2 = - ν e ε φ n j+1 -2φ n j + φ n j-1 (∆x) 2 φ n+1 j -φ n j ∆t + u n j = 0 u 0 j = f j , φ 0 j = g j (10)
• Scheme 2: Semi-implicit (u, φ)-coupling:

           u n+1 j -u n j ∆t -µ u n+1 j+1 -2u n+1 j + u n+1 j-1 (∆x) 2 = - ν e ε φ n j+1 -2φ n j + φ n j-1 (∆x) 2 φ n+1 j -φ n j ∆t + u n+1 j = 0 u 0 j = f j , φ 0 j = g j (11)
• Scheme 3: Implicit (u, φ)-coupling:

           u n+1 j -u n j ∆t -µ u n+1 j+1 -2u n+1 j + u n+1 j-1 (∆x) 2 = - ν e ε φ n+1 j+1 -2φ n+1 j + φ n+1 j-1 (∆x) 2 φ n+1 j -φ n j ∆t + u n+1 j = 0 u 0 j = f j , φ 0 j = g j (12)
In practice Scheme 1 is not used, as the new value of u is naturally available for the resolution of φ. Scheme 2 is rather implemented, with no additional cost. However Scheme 1 is still studied for the purpose of making comparison with litterature conditions.

Linear stability analysis

This section is concerned with the numerical stability of the time-stepping schemes defined above, which means that we expect to detect unbounded oscillations of the numerical solutions as the physical and numerical parameters vary. We look for solutions u n j and φ n j that have the form of Fourier modes sums:

u n j = k∈Z ξ n (k)e 2iπkj∆x , φ n j = k∈Z χ n (k)e 2iπkj∆x
By substituting these expressions into each of the former schemes, we can compute an amplification matrix A(k) defined by:

ξ n (k) χ n (k) = (A(k)) n ξ 0 (k) χ 0 (k) ,
whose spectral radius allows to specify necessary stability conditions, through the Von Neumann condition. The aim of this section is the proof of the following lemma, concerning the schemes 1, 2, and 3:

Lemma 3.1. A necessary condition for the stability of the numerical scheme defined by ( 10) is

∆t ≤ µε ν e ( 13
)
Lemma 3.2. A necessary condition for the stability of the numerical scheme defined by ( 11) is

∆t < (µε + max(µε, √ ν e ε∆x)) ν e (14)
Lemma 3.3. The scheme defined by ( 12) is unconditionally stable.

Proof of Lemma 3.1

Let k ∈ Z. By replacing u n j by ξ n (k)e 2iπkj∆x and φ n j by χ n (k)e 2iπkj∆x in (10), and by removing the k-dependence of χ n and ξ n for the sake of readability, we get

   1 + 4µ∆t (∆x) 2 sin 2 (πk∆x) ξ n+1 = ξ n + 4ν e ∆t ε(∆x) 2 sin 2 (πk∆x) χ n χ n+1 = χ n -∆tξ n (15) We define β(k) = 4ν e ∆t ε(∆x) 2 sin 2 (πk∆x) and δ(k) = 1 + 4µ∆t (∆x) 2 sin 2 (πk∆x) to obtain ξ n+1 χ n+1 =   1 δ(k) β(k) δ(k) -∆t 1   ξ n χ n = A(k) ξ n χ n
The eigenvalues m ± of A(k) are the solutions of

δ(k) 1 + ∆tβ(k) m 2 - δ(k) + 1 1 + ∆tβ(k) m + 1 = 0 As δ(k) = 1 + µε ν e β(k), this equation is equivalent to 1 + µε νe β(k) 1 + ∆tβ(k) m 2 - 2 + µε νe β(k) 1 + ∆tβ(k) m + 1 = 0
whose discriminant is given by:

∆(k) = β(k) (1 + ∆tβ(k)) 2 µε ν e -4∆t µε ν e β(k) -4∆t
Let us note that for k = 0, β(k) = 0. In this case we solve the equation m 2 -2m + 1 = 0, which admits the double root 1. We suppose thereafter that k = 0, which means β(k) > 0 as ν e > 0.

Using that

1 + ∆tβ(k) 1 + µε νe β(k)
stands for the square of the complex root's module when ∆(k) < 0, we distinguish between three cases:

• For ∆t > µε ν e , i.e. 1 + ∆tβ(k) 1 + µε νe β(k)
> 1, there must hold ∆(k) < 0 and we get two complex conjugate eigenvalues of module

1 + ∆tβ(k) 1 + µε νe β(k) > 1.
Therefore the Von Neumann condition is not satisfied.

• For ∆t < µε ν e , i.e.

1 + ∆tβ(k) 1 + µε νe β(k) < 1,
-If ∆(k) < 0, we get two complex conjugate eigenvalues of module

1 + ∆tβ(k) 1 + µε νe β(k) < 1.
-If ∆(k) > 0, we get two real and distinct eigenvalues, which read

m ± = 2 + µε ν e β(k) ± β(k) µε ν e -4∆t µε ν e β(k) -4∆t 2 + 2 µε ν e β(k) Thus 2 + µε ν e β(k) -β(k) µε ν e 2 β(k) 2 + 2 µε ν e β(k) ≤ m ± and m ± ≤ 2 + µε ν e β(k) + β(k) µε ν e 2 β(k) 2 + 2 µε ν e β(k) i.e. 0 < 2 2 + 2 µε ν e β(k) ≤ m ± ≤ 2 + 2 µε ν e β(k) 2 + 2 µε ν e β(k) = 1 Therefore |m ± | ≤ 1.
-Finally if ∆(k) = 0 we get the double root

m = 2 + µε ν e β(k) 2 + 2 µε ν e β(k) Which verifies |m| < 1.
For ∆t < µε ν e the Von Neumann condition is therefore always satisfied.

• Finally for ∆t = µε ν e , i.e. In conclusion, a necessary condition for the numerical stability of Scheme 1 is

1 + ∆tβ(k) 1 + µε νe β(k) = 1, then ∆(k) = β(k) (1 + ∆tβ(k)) 2 -3 µε ν e 2 β(k) -4∆t < 0,
∆t ≤ µε ν e (16)
We point out that for µ = 0, this scheme is then unconditionally unstable.

Proof of Lemma 3.2

By an analogous argument as in the former section, the stability analysis reduces to the computation of the solutions m ± of the following equation:

1 + µε ν e β(k) m 2 -2 + µε ν e -∆t β(k) m + 1 = 0
whose discriminant is given by:

∆(k) = 4∆tβ(k) µε ν e -∆t 2 ν e ε(∆x) 2 sin 2 (πk∆x) -1
We rule out the case k = 0 (double root 1), and we point out that

1 1 + µε ν e β(k) < 1
Let us note w = µε ν e -∆t 2 ν e ε(∆x) 2 . As before we distinguish between three cases:

• If w ≤ 1, i.e. ∆t ∈ µε ν e - √ ε∆x √ ν e ; µε ν e + √ ε∆x √ ν e , then ∆(k) ≤ 0.
For ∆(k) < 0 we get two complex conjugate eigenvalues of module

1 1 + µε ν e β(k) < 1. For ∆(k) = 0 the double root reads m = 2 + µε ν e -∆t β(k) 2 + 2 µε ν e β(k) . So 0 ≤ m ≤ 1.
Thus the Von Neumann condition is satisfied in this case.

• If w > 1 and ∆t < µε ν e -√ ε∆x √ ν e , then there exists some k values such that

∆(k) > 0.
For k such that ∆(k) ≤ 0 we obtain stability by a similar argument as before.

If ∆(k) > 0 instead, the product and the sum of the two eigenvalues are positive, so m ± ≥ 0 > -1. As m -≤ m + , it suffices to prove that m + ≤ 1.

We have

m + ≤ 1 + µε ν e -∆t β(k) 1 + µε ν e β(k) ≤ 1
Hence the Von Neumann condition is satisfied for w > 1 and ∆t < µε ν e -√ ε∆x √ ν e .

• Finally, if w > 1 and ∆t > µε ν e + √ ε∆x √ ν e , as m -≤ m + it remains to prove that m + ≤ 1 and m -≥ -1.

We get readily m

+ ≤ 2 2 + 2 µε ν e β(k) ≤ 1.
Besides,

m -≥ -1 ⇔ 2 + µε ν e -∆t β(k) -β(k) µε ν e -∆t 2 β(k) -4∆t ≥ -2 -2 µε ν e β(k) ⇔ 4 + 3µε ν e -∆t β(k) > β(k) µε ν e -∆t 2 β(k) -4∆t ≥ 0
By taking the square of each side, we get

m -≥ -1 ⇒ µε ν e 2µε ν e -∆t β 2 (k) + 6µε ν e -∆t β(k) + 4 ≥ 0 (17)
Let's recall that β(k) ≥ 0 and study the polynomial in β(k). Its discriminant is:

∆ = 2µε ν e + ∆t
     φ n+1 j -2φ n j + φ n-1 j (∆t) 2 - ν e ε φ n j+1 -2φ n j + φ n j-1 (∆x) 2 = 0 φ 0 j = g j , φ 1 j = g j -∆tf j (19)
This is a discretization of the wave equation

     ∂ 2 φ ∂t 2 - ν e ε ∂ 2 φ ∂x 2 = 0 φ(0, x) = g(x), ∂φ(0, x) ∂t = -f (x) (20)
with an explicit scheme centered in space and time. The time step stability condition for such a scheme is known to reads:

∆t < √ ε∆x √ ν e (21)
We can check that taking µ = 0 in (18) immediately gives (21).

Proof of Lemma 3.3

In order to determine a stability condition for the Scheme 3, we solve the following equation:

1 + µε ν e + ∆t β(k) m 2 -2 + µε ν e β(k) m + 1 = 0
whose discriminant is given by:

∆(k) = β(k) µε ν e 2 β(k) -4∆t
As before, we rule out the case k = 0 (double root 1), and we point out that

1 1 + µε ν e + ∆t β(k) < 1
For ∆(k) < 0 we get two complex conjugate eigenvalues of module

1 1 + µε ν e + ∆t β(k) < 1
When ∆(k) > 0, the two real eigenvalues read

m ± = 2 + µε ν e β(k) ± β(k) µε ν e 2 β(k) -4∆t 2 + 2 µε ν e + ∆t β(k) So 2 + µε ν e β(k) -β(k) µε ν e 2 β(k) 2 + 2 µε ν e + ∆t β(k) ≤ m ± ≤ 2 + µε ν e β(k) + β(k) µε ν e 2 β(k) 2 + 2 µε ν e + ∆t β(k) Thus 0 ≤ 2 2 + 2 µε ν e + ∆t β(k) ≤ m ± ≤ 2 + 2 µε ν e β(k) 2 + 2 µε ν e + ∆t β(k) ≤ 1 Therefore |m ± | ≤ 1. Finally, if ∆(k) = 0 the double root is m = 2 + µε ν e β(k) 2 + 2 µε ν e + ∆t β(k) Thus |m| ≤ 1.
In conclusion, the Von Neumann necessary condition is always satisfied for Scheme 3.

Comments

According to the Lemma 3.1, 3.2, and 3.3, we claim that Scheme 1 and 2 are conditionally stable, whereas Scheme 3 is unconditionally stable. Despite the conditions ( 16) and ( 18) are only necessary ones, the relation between viscous and elastic effects clearly appears. Indeed, the strongest time step restriction occurs for a small fluid viscosity and a large elastic force. In this respect these conditions are in good accordance with the condition computed in [START_REF] Vigneaux | Méthodes Level Set pour des problèmes d'interface en microfluidique[END_REF] with a mathematical analysis of the continuous model. More precisely, the parameter ε in ( 16) and (18) stands for the size of an interface's neighbourhood on which the elastic force is applied. In practice, this parameter is taken in the order of ∆x. By taking ε = ∆x in ( 16) and ( 18) we obtain

∆t ≤ µ∆x ν e (22)
for Scheme 1, and

∆t < µ∆x + max(µ, √ ν e ∆x)∆x ν e (23)
for Scheme 2.

We recognize for Scheme 1 the condition of [START_REF] Vigneaux | Méthodes Level Set pour des problèmes d'interface en microfluidique[END_REF]. Moreover, the semi-implicit treatment of coupling in Scheme 2 allows us to relax easily the former condition. Indeed, for ∆x and µ fixed, beyond a certain value of ν e , we reach a 1/ √ ν e -dependence for the time step, instead of the 1/ν e -dependence related to the explicit coupling. Therefore this condition is by far less restrictive when ν e is large. As the scheme studied in [START_REF] Vigneaux | Méthodes Level Set pour des problèmes d'interface en microfluidique[END_REF] is the same as Scheme 2, we may explain this relaxation on the condition by the fact that the mathematical study performed in [START_REF] Vigneaux | Méthodes Level Set pour des problèmes d'interface en microfluidique[END_REF] did not benefit fully from the semi-implicit coupling between u and φ. We can confirm this argument by the fact that the condition of [START_REF] Vigneaux | Méthodes Level Set pour des problèmes d'interface en microfluidique[END_REF] matches (22) instead.

Moreover, when the fluid viscosity is small, (23) reads approximately

∆t < ∆x 3 ν e (24)
We then retrieve a condition similar to that computed in [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF] by an heuristic method, with a dependence in ∆x In conclusion the stability condition relative to Scheme 2 is unifying the conditions of [START_REF] Vigneaux | Méthodes Level Set pour des problèmes d'interface en microfluidique[END_REF], valid when viscous effects are large and the one computed in [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF] when inertial effect are dominating.

Numerical tests

In this section we numerically study the stability of the following scheme:

                         u n+1 -u n ∆t + u n .∇u n -µ∆u n+1 + ∇p n+1 -F (φ n ) = 0 F (φ n ) = E ′′ (|∇φ n |) ∇ 2 φ n ∇φ n |∇φ n | ∇ × φ n |∇φ n | ∇ × φ n |∇φ n | -E ′ (|∇φ n |)κ(φ n )∇φ n ] 1 ε ζ φ n ε|∇φ n | div u n+1 = 0 φ n+1 -φ n ∆t + u n+1 ∇φ n = 0
This is a 2D non linear version of Scheme 2. For this model we enforce the following condition on time step:

∆t ≤ ∆t A , ∆t A = C (µ + max(µ, √ ν e ∆x))∆x ν e (25)
with C a positive constant determined numerically.

We validate numerically the condition (23) on a simple test case currently used in the litterature ( [START_REF] Lee | Immersed interface methods for incompressible flow with moving interfaces[END_REF], [START_REF] Boffi | Numerical stability of the finite element immersed boundary method[END_REF], [START_REF] Newren | Enhancing the Immersed Boundary Method: stability, volume conservation, and implicit solvers[END_REF]): an elastic membrane, which has initially an elliptic shape, relaxes toward a circular equilibrium state. The membrane is immerged in an incompressible fluid so that the interior area it delimits remains constant. The magnitude and frequency of oscillations around the equilibrium form depends on the fluid viscosity and the membrane's stiffness.

The initial membrane is an ellipse with major and minor axes ã = 0.65 and b = 0.575 respectively. Its stretching, constant along the curve, is e = 1.2526. This configuration corresponds to an equilibrium circle of radius R ≃ 0.6113. The initial and equilibrium states are then quite close as the gap between the final radius and the initial axes is up to 6%. The initial velocity and pressure are set to zero. The computational domain is [0, 2] × [0, 2], and we set ρ = 1 and µ = 1 throughout the domain. We define M = 2/∆x.

Our aim is to show that our stability condition is sharper than the two other conditions on the classical test described above when we vary coefficients. Thus in a first test, we fix ∆x and test the dependence in ν e of the condition (23), for ν e from 1 to 10 6 , with ∆x = 0.03125 (M = 64). We control the time evolution of the ellipse axes. Each simulation is displayed for t ∈ [0, 200/ √ ν e ]. The time steps ∆t A used in the simulations, obtained from (23) are presented in table 1, in comparison with the maximum time steps recommended by the conditions of [START_REF] Vigneaux | Méthodes Level Set pour des problèmes d'interface en microfluidique[END_REF] and [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF]. The same constant C is used for ∆t A and ∆t V even if we could adjust it. However, the goal of the comparison is mostly to show the discrepancy between the time steps ∆t A and ∆t V as ν e grows up. The maximum time step ∆t V , valid for small ν e (with respect to viscosity), becomes far too much restrictive for high ν e . Table 1: The largest time steps allowed for stability, for µ = 1, M = 64 and ν e from 1 to 10 6 : ∆t A computed with the condition (23) and used in the simulations, ∆t V and ∆t B determined with the conditions of [START_REF] Vigneaux | Méthodes Level Set pour des problèmes d'interface en microfluidique[END_REF] and [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF].

From the pictures 1 and 2 we can deduce that the condition (23), less restrictive than the condition of [START_REF] Vigneaux | Méthodes Level Set pour des problèmes d'interface en microfluidique[END_REF] as ν e grows up, enables to ensure stability for a certain range of stiffness (ν e ≤ 10 4 ). Note that for ν e = 10 6 , the equilibrium state is not reached, but for this kind of stiffness, strong non-linear instabilities are expected, and thus the conclusions of the linear analysis are no longer significant (this is in contrast with an unstable case with a time-step higher that ∆t A for which the simulation blows up in a few steps). Moreover, the oscillations size decreases quite rapidly, and becomes smaller that the grid size ∆x = 0.03125, which may be another cause of instability. A last clue that the instability for ν e = 10 6 has nothing to do with linear instability is that it persists when the time step is divided by e.g. 10. This suggests that an implicit coupling may fail to guarantee stability, in spite of an increased computational cost. Indeed, computations are performed with a massless membrane, which means that the fluid density is set everywhere. Thus we are typically in the case of a large added mass effect, with strong interactions between pressure field variations and the membrane. It is shown in [START_REF] Borazjani | Curvilinear immersed boundary method for simulating fluid structure interaction with complex 3D rigid bodies[END_REF] that in this case the convergence of the iterative method in order to solve each time step may be hard to be reached, unless an under-relaxation scheme with a problem dependent coefficient is used.

Next we study the ∆x-dependence of the condition (23), with ν e = 1 fixed. As before we present in table 2 the time steps ∆t A , ∆t V , and ∆t B , for M = 64, 128, 256. This is to advocate that the ∆x Table 2: The largest time steps allowed for stability, for µ = 1, ν e = 1 and M from 64 to 256: ∆t A computed with the condition (23) and used in the simulations, ∆t V and ∆t B determined with the conditions of [START_REF] Vigneaux | Méthodes Level Set pour des problèmes d'interface en microfluidique[END_REF] and [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF].

The time evolution of the membrane axes displayed on picture 3 allows to confirm the ∆x-dependence of the condition (23).

= 1 .

 1 and we get two complex conjugate eigenvalues of module 1 + ∆tβ(k) 1 + µε νe β(k) Thus the Von Neumann condition is still satisfied.

× 10 - 9 C

 9 ν e ∆t A with C = 0.2 ∆t V with C = 0.2 ∆t B 1 1.25 × 10 -2 6.25 × 10 -3 C × 5.52427 × 10 -3 × 5.52427 × 10 -6

Figure 1 :

 1 Figure 1: The time evolution of the membrane axes for ν e = 1 (left) and ν e = 10 2 (right), µ = 1, M = 64.

Figure 2 :

 2 Figure 2: The time evolution of the membrane axes for ν e = 10 4 (left) and ν e = 10 6 (right), µ = 1, M = 64.

3 2 3 128 6 . 4 256 3 .

 33643 dependence of ∆t B is not sharp when surface tension effects are small.M ∆t A with C = 0.2 ∆t V with C = 0.2 ∆t B 64 1.25 × 10 -2 6.25 × 10 -3 C × 1.10485 × 10 -25 × 10 -3 3.125 × 10 -3 C × 3.90625 × 10 -125 × 10 -3 1.5625 × 10 -3 C × 1.38107 × 10 -4

In conclusion, the dependence of the condition (23) in ∆x and 1/ √ ν e is numerically validated on a non-linear 2D test case, for a certain stiffness range. As far as the ν e -dependence is concerned, the stability of the naturally implemented semi-implicit scheme seems to be less conditioned than it was anounced in [START_REF] Vigneaux | Méthodes Level Set pour des problèmes d'interface en microfluidique[END_REF]. Furthermore, the condition (23) and that proposed in [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF] give rise to time steps of the same magnitude when the viscous effects are small with respect to the surface tension effects.

Summary and conclusions

In this work we studied the numerical stability of three time-stepping schemes for a fluid-elastic membrane coupling model. The linear stability analysis based on a simplified 1D model enabled us to obtain stability conditions for explicit and semi-implicit coupling. According to the flow regime, these conditions can be compared to the conditions proposed either in [START_REF] Vigneaux | Méthodes Level Set pour des problèmes d'interface en microfluidique[END_REF] or in [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF] with different methods. On the one hand, the dependence in 1/ √ ν e of the condition for semiimplicit coupling accounts for a relaxation of the condition of [START_REF] Vigneaux | Méthodes Level Set pour des problèmes d'interface en microfluidique[END_REF], especially for problems with strong membrane stiffness. On the other hand, the dependence in ∆x 1 2 of the same condition allows to compare it to that of [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF]. The numerical tests on a non-linear 2D case validate the analysis results for a certain range of membrane's stiffness. The study of the added mass effect and the implementation of an efficient implicit scheme are under consideration for further work.


Roger S Bagnall

Kai Brodersen

Craige B Champion

Andrew Erskine

Sabine R Huebner

Horapollon

Jean-Luc Fournet

To cite this version: Horapollon is known as the author of a treatise on the HIEROGLYPHS (Hieroglyphica), the only one in all preserved Greco-Latin literature. He seems to belong to a family originating from the Thebaid (Egypt), hypothetically reconstructed with the help of various literary sources and one papyrus: his homonymous grandfather (end fourth/beginning fifth century CE) taught Greek grammar in Alexandria and Constantinople; his father Asklepiades wrote (lost) works on pharaonic history and religion; and his uncle (and father-in-law) Heraiskos was a pagan philosopher. Horapollon himself taught philosophy in the last third of the fifth century in Alexandria and, after being tortured and suffering a reverse of fortune, finally converted to Christianity. So Horapollon was a member of the pagan intelligentsia that attempted to resist Christianity (now the state religion) by exalting the old, secular Hellenism and concomitantly defending the dying pharaonic culture (the last known use of hieroglyphics dates from 394). This biculturality accounts for the Hieroglyphica: while about half of the hieroglyphic signs dealt with by Horapollon are more or less correct (which implies a real, even though superficial and deficient, knowledge of this writing), the exegesis he gives of them, purely allegorical, relies on traditions of late Hellenism.

Though this treatise was of no help to modern scholars (before Champollion) in deciphering hieroglyphics, it nevertheless nourished the artistic repertoire of the Renaissance and the vogue of the emblematic literature.

SEE ALSO: Alexandria (Egypt); Egyptomania, eastern; Egyptomania, western; Hellenism, Byzantine; Paganism, Byzantine.