D Sheeren

N Bastin

A Ouin

S Ladet

G Balent

J.-P Lacombe

Discriminating small wooded elements in rural landscape from aerial photography: a hybrid pixel/object-based analysis approach

While small, fragmented wooded elements do not represent a large surface area in agricultural landscape, their role in the sustainability of ecological processes is recognized widely. Unfortunately, landscape ecology studies suffer from the lack of methods for automatic detection of these elements. We propose a hybrid approach using both aerial photographs and ancillary data of coarser resolution to automatically discriminate small wooded elements. First, a spectral and textural analysis is performed to identify all the planted-tree areas in the digital photograph. Secondly, an object-orientated spatial analysis using the two data sources and including a multi-resolution segmentation is applied to distinguish between large and small woods, copses, hedgerows and scattered trees. The results show the usefulness of the hybrid approach and the prospects for future ecological applications.

Introduction

The importance of green veining in agricultural landscapes has already been demonstrated for soil and biodiversity conservation [START_REF] Grashof-Bokdam | Green veining: landscape determinants of biodiversity in European agricultural landscapes[END_REF]. Green veining, which consists of semi-natural elements such as pond, permanent meadows, field margin, hedgerows and woodlots embedded in crop mosaics, helps to maintain wildlife biodiversity (auxiliary, game) in rural landscapes. Among these semi-natural elements, the identification and characterization of small wooded objects (i.e. copses, hedgerows and scattered trees) are often required in landscape ecology studies: for instance, to predict the response of passerine bird abundance to landscape structure and composition [START_REF] Balent | Modelling bird communities/landscape patterns relationships in a rural area of south-western France[END_REF], or to identify landscape elements and configurations that encourage auxiliary hoverflies [START_REF] Sarthou | Landscape parameters explain the distribution and abundance of Episyrphus balteatus (Diptera, Syrphidae)[END_REF].

Although there is a growing interest in small wooded elements and their role in ecological processes, landscape ecology studies suffer from the lack of operational methods for automatic detection of these elements from aerial photographs. Several investigations have already been carried out regarding the extraction of areas of vegetation or anthropogenic objects (buildings, roads) from remotely sensed data (Baltsavias 2004, Boyd and[START_REF] Boyd | Satellite remote sensing of forest resources: three decades of research development[END_REF]. However, most of the studies related to the identification of wooded elements have been devoted to the classification of large forest patches often with some refinements regarding tree species or boundary *Corresponding author. Email: david.sheeren@ensat.fr delineation (e.g. Peterson et al. 2004, Tuominen and[START_REF] Tuominen | Performance of different spectral and textural aerial photograph features in multi-source forest inventory[END_REF]. The automatic detection of natural, small wooded and thin linear elements in rural landscapes and their differentiation when they are connected but do not underpin the same ecological processes (e.g. adjacent hedges and woods) has drawn less attention [START_REF] Stach | Hedges and tree rows detection with eCognition for the use of the French national forest inventory[END_REF][START_REF] Thornton | Sub-pixel mapping of rural land cover objects from fine spatial resolution satellite sensor imagery using super-resolution pixel-swapping[END_REF][START_REF] Zhang | Automatic extraction of wind erosion obstacles by integration of GIS data, DSM and stereo images[END_REF]. Because of the fine spatial resolution, largescale aerial photographs are preferred to identify such elements by applying computer-aided interpretation techniques. However, this approach is generally time-consuming and, thus, expensive and not recommended when investigating large areas. In addition, the spectral information provided by this kind of data is not always sufficient to obtain classifications of high accuracy. It often requires adding texture into the classification step. Medium and high spatial resolution remotely sensed data (e.g. SPOT and Landsat Thematic Mapper (TM) images) are also widely used to identify wooded areas. These data provide richer spectral information and cover larger areas. However, the relatively coarse spatial resolution (20 m • 20 m) does not enable detection of all small wooded fragments such as copses, hedgerows and scattered trees.

Since the availability of very high spatial resolution satellite images, object-based image analysis (OBIA) has been proposed as an alternative to the pixel-based classification approaches [START_REF] Benz | Multiresolution, object-oriented fuzzy analysis of remote sensing data for GIS-ready information[END_REF]. It involves segmenting images into homogeneous regions and characterizing objects with a set of features related to spectral, spatial and contextual properties. Because these image segments represent objects in the landscape more accurately than the original pixels, and classify them according to their shapes and sizes, this approach is well suited for the extraction of small wooded elements. In addition, previous studies have shown that object-based methods outperform pixel-based methods when applied to aerial photographs (e.g. [START_REF] Cleve | Classification of the wildlandurban interface: a comparison of pixel-and object-based classification using highresolution aerial photography[END_REF]. However, segmentation in objectbased classification is still a crucial task. The user must decide the specific scale at which an image is to be segmented into objects and this scale can vary from one class to another. Scales are often difficult to define and parameter values are highly dependent on the user's experience. In that sense, a hybrid pixel/object-based classification may be less dependent on the results of the segmentation. Such a hybrid approach has already been proposed with success. Among the few recent researches in this area, [START_REF] Wang | Integration of object-based and pixel-based classification for mapping mangroves with IKONOS imagery[END_REF] integrated the pixel and object-based classification for mapping mangroves from IKONOS imagery. [START_REF] Shackelford | A combined fuzzy pixel-based and object-based approach for classification of high-resolution multispectral data over urban areas[END_REF] adopted a hybrid fuzzy approach for urban land-cover classification. Other works presented a similar method for shape-based object recognition: [START_REF] Li | A hybrid object-based/pixelbased classification approach to detect geophysical phenomena[END_REF] for geophysical phenomena detection and Carmichael and Herbert (2004) for computer vision applications. In general, two main conclusions can be drawn from the above studies.

1. A more accurate classification can be achieved with a hybrid approach: a pixelbased classification can be applied for classes with good separability while an object-based analysis can be performed to better discriminate spectrally mixed classes. 2. A more accurate segmentation can also be obtained: the classified pixel-level can serve to optimally segment the individual objects of interest.

In this paper, we propose a hybrid method coupling a pixel-based classification approach to extract woody elements in the landscape using texture from aerial photographs and an object-based image analysis to discriminate among the different wooded element types: large and small woods, copses, hedgerows and scattered trees.

Material

The study area is included in the territory of the Long-Term Ecological Research (LTER-Europe) study site 'Valleys and Hills of Gascogne' (VHG), France (43 13¢ N, 0 52¢ E). It is an agricultural landscape typical of south-western France (agricultural mosaic with wooded patches) with mixed production systems, i.e. livestock farming systems and cropping systems. The total area of the study site is approximately 25 000 ha and includes a meso-scale watershed (12 000 ha) -the Ne `re river.

The major data used for the study were a mosaic of orthorectified digital aerial true colour photographs acquired in 2002 (named BDORTHO Ò ) and produced by the French mapping agency (IGN). This set of geo-referenced photographs covers approximately 31 km 2 and has a spatial resolution of 0.5 m • 0.5 m. Infrared (IR) orthophotos which are more appropriate for extracting wooded elements also exist for this study site. However, these data have not been used, the most recent date being 1996 which is too remote to combine with the present ecological data. In addition, the structure of the agricultural land has since been modified and has led to the disappearance of some hedges.

Additional data were also used to discriminate wooded areas during the objectorientated analysis, especially to disconnect hedgerows from wood patches. These data are a binary wood/non-wood image derived from an existing supervised landcover classification [START_REF] Barreau | Dynamique spatio-temporelle d'Episyrphus balteatus dans la fore ˆt fragmente ´e des coteaux de Gascogne[END_REF]. This classification was based on four SPOT-4 multi-spectral (MS) multi-temporal images acquired in 2001 (January, April, July and October) with a spatial resolution of 20 m • 20 m. In a general way, any source of data of the same period providing a wood/non-wood map with coarser resolution (i.e. .10 m • 10 m) can be used as ancillary data in the method (e.g. the CORINE Land Cover map or the Pan-European Forest/non-Forest map proposed by the EU-Joint Research Centre [START_REF] Pekkarinen | Pan-European Forest/Non-Forest mapping with Landsat ETM+ and CORINE Land Cover 2000 data[END_REF]).

Method

The hybrid method we propose has two main steps: (1) a spectral and textural pixelbased analysis to get all the wooded elements of the study area and (2) a spatial objectbased analysis to delineate the different wooded element types: small and large woods, copses, scattered trees and linear elements (figure 1). Data pre-processing was performed to ensure accurate geometrical correspondence between the two sources of data: the true-colour photograph and the binary wood/non-wood mask extracted from the available land-cover classification. In addition, the digital photograph was resampled to a pixel size of 1 m • 1 m for optimal running time computation during textural analysis and segmentation. The pixel-based image analysis was conducted using the resampled digital photograph mosaic.

Pixel-based image analysis: extraction of wooded elements

3.1.1 Textural analysis. Texture is the expression of local variations in grey-levels in the pixels of the digital image that result from the particular organization of objects in the landscape [START_REF] Haralick | Texture features for image classification[END_REF]. Several studies demonstrated interest in incorporating such spatial information into classification procedures to improve classification accuracy. The addition of texture has become quite common for land-cover mapping [START_REF] Marceau | Evaluation of the grey-level co-occurrence matrix-method for land cover classification using SPOT imagery[END_REF][START_REF] Puissant | The utility of texture analysis to improve perpixel classification for high to very high spatial resolution imagery[END_REF]. It has been also widely used in forestry and ecological mapping applications [START_REF] St Onge | Estimating forest stand structure from high resolution imagery using the directional variogram[END_REF][START_REF] Franklin | Incorporating texture into classification of forest species composition from airborne multispectral images[END_REF][START_REF] Coburn | A multiscale analysis procedure for improved forest stand classification[END_REF], especially when digital aerial photographs are required [START_REF] Wulder | Aerial image texture information in the estimation of northern deciduous and mixed wood leaf area index (LAI)[END_REF][START_REF] Caridade | The use of texture for image classification of black & white air photographs[END_REF]. Generally, the resolution of aerial photographs is finer than or similar to the size of the trees in the image. Thus, local variations in the pixels can be related to forest structure variables [START_REF] Tuominen | Performance of different spectral and textural aerial photograph features in multi-source forest inventory[END_REF]Pekkarinen 2005, Kayitakire et al. 2006).

In this study, the textural information is not used for estimating wood attributes (e.g. species, height and density). The texture is used mainly to extract all the wooded elements from the orthophotos more efficiently and limit confusion with other classes that have similar pixel values (e.g. grasslands and crop fields). Since the further supervised classification step can be performed without this information, the textural analysis -which is time consuming -can be viewed as an optional step in the method. However, as illustrated in the results, adding this information contributes to considerable improvement in the classification accuracy.

There are several approaches to texture processing. The most widely used in remote sensing is based on the grey-level co-occurrence matrix (GLCM) proposed by [START_REF] Haralick | Texture features for image classification[END_REF]. This matrix computes the relationships of pixel values from which various second-order statistical texture measures can be derived [START_REF] Haralick | Statistical and structural approaches to texture[END_REF]. Several parameters must be defined for computation of the GLCM: the size of the moving window, and offset distance and direction.

In our case, height texture measures were computed from the green band of the aerial photograph: dissimilarity, contrast, homogeneity, variance, entropy, correlation, angular second moment and mean. This set of measures was chosen with the intent of not selecting arbitrarily just a few of them. The 8-bit image (256 grey levels) was rescaled into 32 grey levels to ensure statistical validity of the measures [START_REF] Clausi | An analysis of co-occurrence texture statistics as a function of grey-level quantization[END_REF]) and to reduce the computation time of the GLCM (which is proportional to 2 n where n is the number of bits). Four main directions (0 , 45 , 90 , 135 ) and three window sizes (3 • 3, 7 • 7, 11 • 11) were used with a constant offset distance of one pixel. The combination of all GLCM parameter values produced 96 texture images.

Figure 1. Flowchart of the method proposed to discriminate small wooded elements in rural landscape. A pixel-based image analysis is performed to extract wooded objects from truecolour photography. An object-based analysis is applied in the next step using the result of the previous step and an additional wood/non-wood map with coarser resolution, to enhance the thematic resolution of the pixel-based classification and disconnect objects related to different ecological processes.

Finally, to reduce the number of the texture images computed and decorrelate them without decreasing the main textural information, standardized principal component analyses (PCA) based on correlation matrices were performed following two different strategies (figure 2). The strategies were tested to compare their performance and to identify a possible optimal window size for wooded objects with the understanding that the success of the classification procedure using textural measures depended largely on the selected window size [START_REF] Puissant | The utility of texture analysis to improve perpixel classification for high to very high spatial resolution imagery[END_REF]). In the first strategy (figure 2(a)), a PCA was applied on all the texture variables including the different window sizes and directions. In the second strategy (figure 2(b)), a PCA was computed for each set of texture variables related to only one window size including each direction. With each strategy, more than 90% of the textural information was conserved within the principal components retained for the classification.

3.1.2 Classification and post-processing. Several classifications were performed following the different strategies (see C 1 , C 2 , C 3 , C 4 in figure 2). The principal components were added to the R-G-B bands of the digital photograph to extract all the wooded areas and obtain a binary mask (i.e. the planted/non-planted-tree mask mentioned in the next section). The standard maximum likelihood classifier (MLC) was applied using training site data derived from visual interpretation. In the next post-classification step, the resulting mask was filtered to reduce the salt-and-pepper effect in the image and to join some disconnected planted-tree elements within the woods. First, a major analysis with a window size of 9 • 9 was performed on pixels classified as planted-tree areas. Then, the same filter was applied successively on nonplanted-tree areas in an iterative way until all the wood patches were filled in. This selective post-processing enabled some classification artefacts to be corrected without Figure 2. The classification strategies tested to assess the usefulness of adding texture. Each cell of the matrix represents one texture image related to one texture measure (DISS, dissimilarity; CON, contrast; HOM, homogeneity; ASM, angular second moment; ENT, entropy; COR, correlation; VAR, variance; MEAN), one window size and one direction. In the first strategy (a), a principal component analysis (PCA) was computed using all the texture variables including all the window sizes and directions. In the second strategy (b), a PCA was applied on each set of texture variables related to one window size and including all directions. substantial alteration of the borders of the wooded elements. Finally, classification C 1 , which proved to be the best, was retained for the following object-orientated image analysis.

Object-based image analysis: discrimination between wooded element types

A spatial analysis based on the object-orientated approach was performed to get more details in the planted-tree areas class and to disconnect some merged elements that were semantically different and related to different ecosystems (in particular, the woody fragments and the hedgerow network). This analysis was carried out following the OBIA methodology composed of three steps: segmentation, object characterization, classification. These steps are described in detail in the following sections.

3.2.1 Multi-resolution segmentation. Image objects were created from the binary mask obtained in the previous step using the multi-resolution segmentation method proposed by Definiens Software (eCognition). The segmentation algorithm introduced by Baatz and Scha ¨pe ( 2000) is a bottom-up region-growing technique which merges stepwise adjacent pixels according to their similarity. Two pixels are merged to form larger segments if the resulting segments do not exceed a user-defined heterogeneity threshold. The definition of heterogeneity includes both spectral and spatial criteria and is fixed by the 'scale parameter'. Distinct weights can be assigned to the data layers used in the segmentation process.

In this study, segmentation was conducted at two levels using different scale parameters to construct a hierarchical image object network. The first level (level 1) was created with a scale parameter fixed at 200 using only the spectral heterogeneity criterion related to the binary mask composed of the planted/non-planted-tree classes (weight = 1). The R-G-B bands of the digital photograph were not used at this level (weight = 0). In this way, all the pixels classified in the binary mask were aggregated and transformed into image objects with boundaries strictly identical to the classes at the pixel level. In the second step, a lower level (level 2) was created with a scale parameter of 50 using both spectral and spatial heterogeneity criteria. This level was processed from the three visible bands of the aerial photograph (weight = 1) but without using the binary mask of planted-tree areas (weight = 0). The sub-objects created at this level with finer segmentation enabled disconnected elements to be merged at the higher level (e.g. adjacent woods and hedges).

It is important to note that the method is not very sensitive to the scale parameter values defined for the multi-resolution segmentation. At segmentation level 1, the scale value must be sufficiently high to avoid an over-segmentation of the pixel-based classification. This condition does not cause any difficulty since the segmentation is performed on the binary mask of the planted/non-planted area classes. On the contrary, at segmentation level 2, the scale value must be relatively low to over-segment the aerial photograph and disconnect elements belonging to different classes.

Class definition and image object characterization.

The image objects were characterized with different features according to the classes we wanted to extract. The class hierarchy defined for this study is presented in figure 3. For level 2 classification, we defined two classes of woody areas: classes 'A' and 'B'. The first included the wooded elements classified as planted-tree areas in both the data sources (i.e. the planted-tree areas of the level 1 segmentation derived from the aerial photograph and the wooded/non-wooded areas derived from the existing land-cover classification of the 20 m • 20 m resolution SPOT images). Thus, this class 'A' contained mainly large planted-tree areas but not the small fragments of woods and hedges that did not appear in the data with coarser resolution. These small objects were represented in class 'B', which included the planted-tree areas that only existed in the aerial photograph (i.e. in segmentation level 1). The definitions adopted for the child classes at level 2 were adapted from the French National Forest Inventory (IFN) typology. A woody area is classified as a 'large wood' if it covers more than 4 ha. The 'small wood' class includes areas of between 50 ares and 4 ha. The 'copse' class includes areas of between 5 and 50 ares. 'Scattered trees' are clumps of less than 5 ares and 'natural linear elements' are tree rows with a width of less than 25 m and a length of more than 25 m. In contrast to the IFN definitions, we did not take into account the average width at the top of the trees or the species composition in our definitions.

According to this typology, several spatial and relational features were computed to characterize the objects and to prepare the rule-based classification. The features selected were: area, length, width, compactness, elliptical fit and relation to superobjects. The wood/non-wood mask with coarser resolution derived from the landcover classification was imported as a raster thematic layer. The affiliation of the image objects to the thematic class of this layer was defined through a thematic attribute.

Rule-based classification.

A rule-based classification composed of several steps concluded the object-orientated analysis. First, the image objects at segmentation level 1 were directly classified as planted-tree and non-planted-tree areas using the mask derived from the pixel-based analysis. Next, the image objects at the lower segmentation level (level 2) were processed to distinguish the woody areas 'A' and 'B' within the planted-tree areas. Initially, woody areas labelled 'A' were identified using both the result of the previous classification (through the 'relational to super-object' feature) and the wood/non-wood thematic layer of the land-cover classification (through a corresponding 'thematic' attribute). Woody areas labelled 'B' were extracted as a result. After this classification, the level 2 segmentation (in which the objects were over-segmented) was modified by merging all adjacent objects of the same class (classification-based fusion). This merging was undertaken to take into account the area criteria in the object recognition. The classification of the woody Figure 3. Image object class hierarchy based on the IFN typology. Woody areas 'A' are plantedtree areas in both the aerial photograph and the wood/non-wood mask with coarser resolution derived from an existing land-cover classification. Woody areas 'B' include the planted-tree areas that exist only in the photograph. areas 'A' and 'B' was then broken down into specialized classes (large wood, small wood, copse, scattered trees and linear elements) according to rules that had been defined previously. This classification procedure was followed by a post-processing step to improve the results of the classification. Except for hedges, adjacent objects of the same or different classes (i.e. from different woody areas 'A' or 'B') were merged by keeping as class value the value of the largest woody object. For instance, if an object classified as copse 'B' was adjacent to an object classified as small wood 'A', the resulting object was integrated in the small wood class 'A' after merging. If the area of the newly created object exceeded the maximum value associated with the definition of the resulting class, the object was integrated in the next higher-level woody class. In the previous example, the new object would be integrated in the large wood class 'A' if it was bigger than 4 ha. At the end of the process, the objects in class 'A' and 'B' were no longer distinguished (figure 4).

Validation

A visual interpretation of the digital photography was performed to obtain data for validation of the pixel and object classifications. More than 1% (i.e. 330 199 pixels) of the classified area was randomly selected to assess the accuracy of the pixel-based classification and 107 segments were retained in the case of the refined object-based classification. Accuracy assessment was based on the confusion matrix and associated indicators: overall accuracy and Kappa coefficient.

Classification results and discussion

A statistical summary of the results is presented for both the pixel-based and objectbased classifications (table 1). In a general way, the gain in thematic resolution with the object-based classification led to a loss in overall accuracy (from 92.6% for the best classification in the initial pixel-based analysis to 83.5% at the end of the object-based classification), due to the creation of some mixed objects during the segmentation. However, because the pixel-based and object-based analyses are sequential, the results of the two approaches should not be compared. In the case of the pixel-based analysis, the accuracy assessment is related to the extraction of the planted-tree/non-plantedtree areas from the aerial photographs, with regard to textural information. For the object-based analysis, because the limits of the planted-tree area class have been definitively fixed at the previous step, only confusion errors involved in the discrimination between the wooded element types within the planted-tree area class are computed. The object-based classification results do not indicate confusion between the wooded element classes and the non-planted-tree areas.

In terms of proportions, the 'large wood' class represents 10.90% (342.53 ha) relative to the entire study area versus 3.27% (102.86 ha) for the 'small wood' class, 2.61% (82.21 ha) for the 'copse' class, 1.14% (35.90 ha) for the 'scattered trees' class, 3.44% (108.20 ha) for the 'linear elements' class and 78.63% (2471.08 ha) for nonplanted-tree areas.

Concerning the pixel-based analysis, the addition of textural information always represented an important improvement in classification accuracy compared with the results obtained with the R-G-B classification alone. Without using the textural information, the overall accuracy was 79.5% with a k coefficient of 0.60 (see C 0 in table 1(a)). There was considerable confusion between grasslands and cropped fields. Thus, even if the textural analysis is an optional step in the method, it is recommended to use this information to increase the accuracy degree of classification. Moreover, the pixel-based classification strategy that provided the highest accuracy included all the texture images in the PCA (i.e. C 1 derived from nine principal components and the R-G-B bands) but the results were very similar to C 3 and C 4 (the difference was less than 1%). It seems difficult to determine the most appropriate classification of these three cases. Concerning window size, the intermediate 7 • 7 dimension gave the most accurate results (C 3 ). The smaller dimension (3 • 3) was more suitable for the detection of small objects (scattered trees and thin linear elements) but shadows and planted-tree areas were mixed. The larger dimension (11 • 11) provided a classification that respected the object boundaries better but there were more omissions of small elements. Since the objects to be extracted are of different sizes, an intermediate window size may be optimal (C 3 ) or alternatively, a strategy that incorporates the three window sizes in the PCA (C 1 ), but this is more time-consuming. Finally, although the results of the classifications were not the same, overall accuracy was relatively high for each and provided a binary mask of planted/non-planted-tree areas of good quality. Concerning the object-based analysis, the results derived from the confusion matrix are given in table 1(b). For large and small woods, the classification was very accurate since the only criterion that distinguished the objects was area. However, for copses, linear elements and scattered trees, accuracy was moderate. There was confusion between classes due to segmentation errors. In some cases, hedgerows were difficult to extract because the linear character could not always be verified, especially when the hedges were interconnected. This created only one image object characterized by a complex shape which did not satisfy the previously defined classification rule. The shape property values get skewed. However, this issue is mitigated by ecological requirements. From the fauna point of view, hedges may be connected and form a continuous network throughout the landscape. The main challenge lies in distinguishing woods and hedges which do not underpin the same ecological processes (e.g. breeding versus mobility). Some adjacent hedges and wood patches are connected to form only one object which can introduce a confusion error. However, the use of the wood/non-wood mask with coarser resolution deriving from the existing land-cover classification means this aggregation effect can be limited. This is the main interest of using ancillary data in this hybrid method. At the lower segmentation level, the largest woody areas (class 'A') were classified using the two data sources. The smallest objects (including hedges) that did not appear in the 20 m land-cover classification because they were too small were identified in the digital photograph. Since the photograph is over-segmented at this segmentation level 2, adjacent hedges are not included in the class 'A' woody areas that border them. They can therefore be individualized (figure 5). The remaining aggregation errors mainly appeared between copses, scattered trees and hedges that were only extracted from the aerial photograph but, as mentioned previously, these errors are less critical because the objects concerned are related to the same ecosystems. In all cases, even when two sources can be used to classify the objects, all the object boundaries will be derived from the digital photograph (i.e. the data source with the greatest geometrical accuracy).

Conclusion

The aim of this study was to assess the potential usefulness of combining pixel-based and object-based image analysis using multi-source data to extract wood patches and linear elements in agricultural landscapes. Compared to the traditional per-pixel classification approach, the hybrid method enables the use of spatial and relational features during the classification process. This helps to discriminate between the different wooded element types based on their shape and size (large and small woods, copses, hedgerows and scattered trees). Compared to object-based image analysis, this approach contributes to classification results that are less dependent on the results of segmentation. Considering classification accuracy, the incorporation of texture measures in the pixel-based analysis proved to be the most reliable way to detect plantedtree areas from aerial photographs. The object-orientated classification could probably be improved by incorporating an intermediate segmentation level to reduce errors due to confusion between copses, scattered trees and linear elements. The assessment of the accuracy of this method could also be extended since the confusion matrix gives only information on the thematic accuracy but not on the geometrical accuracy of the image objects (i.e. the quality of segmentation). The resulting maps will now be used to calculate the fragmentation of wooded patches and connectivity indices. The morphological operators proposed by [START_REF] Vogt | Mapping spatial patterns with morphological image processing[END_REF] to define forest spatial patterns are now being applied for this purpose.

Figure 4 .

 4 Figure 4. Excerpt of the pixel and object-based classifications.

Figure 5 .

 5 Figure 5. Disconnecting adjacent hedges and woods using a layer of woods with coarser resolution.

Table 1 .

 1 Summary statistics of the classification results.

	(a) Accuracy of the pixel-based analysis		(b) Accuracy of the object-based analysis (%)
		C 0	C 1	C 2	C 3	C 4		LW SW	C	L	ST
	OA (%): 79.50 92.61 86.90 92.13 91.72	LW	100 0.00	0.00 0.00 0.00
	Kappa	0.60 0.85 0.74 0.84 0.83	SW	0.00 91.30	0.00 3.03 0.00
							C	0.00 0.00 89.47 24.24 9.09
							L	0.00 8.70 10.53 63.64 9.09
							ST	0.00 0.00	0.00 9.09 81.82
							OA (%)	81.3		Kappa 0.76
	C 0 , without textural information; C 1-4 , with textural information; OA, overall accuracy; LW,
	large wood; SW, small wood; C, copse; L, linear element; ST, scattered trees.


Hamdi Raïssi
email: hamdi.raissi@insa-rennes.fr.

COMPARISON OF PROCEDURES FOR FITTING THE AUTOREGRESSIVE ORDER OF A VECTOR ERROR CORRECTION MODEL

Keywords: Vector error correction model, Model selection, Cointegration, Information criteria, Autocorrelation tests

This paper investigates the lag length selection problem of a Vector Error Correction Model (VECM) by using a convergent information criterion and tools based on the Box-Pierce methodology recently proposed in the literature. The performances of these approaches for selecting of the optimal lag length are compared via Monte Carlo experiments. The effects of misspecified deterministic trend or cointegrating rank on the lag length selection is studied. Noting that processes often exhibit nonlinearities, the cases of iid and conditionally heteroscedastic errors will be considered. Strategies which can avoid misleading situations are proposed.

Introduction

The Vector Error Correction Models (VECM) are often used in the statistical analysis of nonstationary variables since they allow to describe several features. In particular the cointegration analysis of equilibrium relationships between variables is much considered in theoretical research. The dominant test for determining the number of equilibrium relationships, the cointegrating rank, is the Likelihood Ratio (LR) test developed by [START_REF] Johansen | Statistical analysis of cointegration vectors[END_REF][START_REF] Johansen | Estimation and hypothesis testing of cointegration vectors in gaussian vector autoregressive models[END_REF]. Noting that this test depend on the specification of the deterministic part of the VECM, [START_REF] Johansen | The role of the constant and linear terms in cointegration analysis of non-stationary variables[END_REF] proposed LR tests based on his likelihood procedure for testing restrictions on the deterministic parameters. However it is well known that the likelihood inference proposed by Johansen for the analysis of the long run relationships and the deterministic terms strongly depend on the choice of the lag length. Indeed if the short run dynamics are over specified this can entail a loss of efficiency in our multivariate framework since a large number of parameters are introduced in this case. Some authors found that the LR test for the cointegrating may suffer from a substantial loss of power in such a case (see e.g. [START_REF] Boswijk | Dynamic specification and cointegration[END_REF]). If the short run dynamics are under specified the residuals become autocorrelated and the asymptotic theory underlying the Johansen's procedure breaks down [START_REF] Johansen | Likelihood-Based Inference in Cointegrated Vector Autoregressive Models[END_REF], Theorem B.13 p 251). Hence the LR test for the cointegrating rank is not valid even asymptotically in this case. One can expect the similar effects on the LR tests for deterministic parameter restrictions when the autoregressive order is misspecified. Therefore it clearly appears that the choice of the lag length of a VECM is crucial.

When a cointegration analysis is conducted the autoregressive order is usually chosen by considering an information criterion. For instance [START_REF] Cavaliere | Testing for co-integration in vector autoregressions with non-stationary volatility[END_REF] studied long run relationships between interest rates and used the Bayesian Information Criterion (BIC) to fit an autoregressive order to the data. Several papers have been devoted to the study of the performances of information criteria to find the optimal lag length of a VECM. Reference can be made to Ho and Sorensen (1996), [START_REF] Gonzalo | Specification via selection in vector error correction models[END_REF] or [START_REF] Hacker | Optimal lag-length choice in stable and unstable VAR models under situations of homoscedasticity and ARCH[END_REF]. These simulation investigations highlight the difficulty of choosing the autoregressive order of a VECM although the trend parameters are not taken into account in these studies. Numerous information criteria are available for the choice of the autoregressive order of VECM, as for instance the asymptotically efficient Akaike Information Criterion (AIC) or the consistent BIC. A consistent information criterion is such that the selected lag length converge to the true autoregressive order while efficiency relies on the optimal prediction error. Since in our case the main objective is to find the optimal lag length for the VECM, we restrict our attention to the commonly used BIC. Note that considering other consistent information criteria would lead to similar general conclusions.

On the other hand tests based on the autocorrelations of the residuals for checking the adequacy of the lag length in the framework of cointegrated variables have been recently proposed in the literature by [START_REF] Duchesne | Testing for serial correlation of unknown form in cointegrated time series models[END_REF], Brüggemann, [START_REF] Brüggeman | Residual autocorrelation testing for error correction models[END_REF] and [START_REF] Raïssi | Autocorrelation based tests for vector error correction models with uncorrelated but nonindependent errors[END_REF]. Some of these tests are implemented in the software JMulTi. In this paper we compare the lag length selection properties of an consistent information criterion and of the Box-Pierce methodology with possibly misspecified cointegrating rank or the deterministic terms. We focus on the portmanteau tests since the Box-Pierce methodology is routinely used for checking the adequacy of fitted models to time series. Using the properties of each of the methods for specifying the short run dynamics of a VECM, strategies which can reduce misleading choices for the lag length are presented. These strategies are proposed taking into account of possibly misspecified long run relations or deterministic terms, and also the possible presence of nonlinearities.

The remainder of the paper is organized as follows. Section 2 outlines the framework of our study and introduces the deterministic parameters restrictions. The estimation procedure for the VECM is described in Section 3. The tools used for model selection in the Monte Carlo experiments are also presented. In Section 4 the framework of our simulation experiments is given and the results are discussed. Strategies exploiting the properties of the Box-Pierce methodology and of information criteria lag length selection procedures are presented.

The following general notation is used throughout the paper. Considering a d × r dimensional matrix A, we define the orthogonal complement A ⊥ , which is a full column rank d × (dr) matrix such that A ′ A ⊥ = 0. For a given random variable a t we define a t q = (E a t q ) 1/q , where . denotes the Euclidean norm. The trace of a square matrix is denoted by Tr and the determinant by | . |.

Characterization of the model

Let us consider the following VECM

∆x t = µ 0 + µ 1 t + αβ ′ x t-1 + p0-1 i=1 Γ i ∆x t-i + ǫ t , (2.1)
where the process (x t ) is d-dimensional and ∆x t := x tx t-1 . The parameters α and β are of full column rank and of dimension d×r 0 . The Γ i 's are d×d short run parameters, and when p 0 = 1 the sum in (2.1) vanishes. The error process (ǫ t ) is commonly assumed iid Gaussian with positive definite covariance matrix Σ ǫ and such that Eǫ t = 0. However the standard assumption of iid Gaussian errors cannot take into account of nonlinear dynamics which often arise in practice. This framework is then considered to be not realistic in many situations. For instance [START_REF] Trenkler | The polish exchange rate system: A unit root and cointegration analysis[END_REF], [START_REF] Koutmos | Asymmetric volatility transmission in international stock markets[END_REF] or [START_REF] Hiemstra | Testing for linear and nonlinear Granger causality in the stock price-volume relation[END_REF] studied financial variables and found a strong evidence of nonlinear dependence in the data. Numerous models in the literature produce processes with conditional heteroscedasticity as for instance hidden Markov models (see e.g. [START_REF] Amendola | Concepts and tools for nonlinear time series modelling[END_REF]) or MGARCH models (see e.g. [START_REF] Bauwens | Multivariate GARCH models: A survey[END_REF]). Taking into account the dependence of the errors is important for testing the adequacy of linear models as pointed out by [START_REF] Francq | Diagnostic checking in ARMA models with uncorrelated errors[END_REF] or [START_REF] Francq | Multivariate portmanteau test for autoregressive models with uncorrelated but nonindependent errors[END_REF]. In addition the study of VECM in nonstandard situations has attracted much attention in the recent years. Therefore we also consider in our simulations the case of dependent but uncorrelated errors which are such that the results of the Granger representation theorem hold (see equation (2.3) below) as in [START_REF] Rahbek | ARCH innovations and their impact on cointegration rank testing[END_REF], [START_REF] Raïssi | Testing the cointegrating rank with uncorrelated but dependent errors[END_REF] or [START_REF] Hacker | Optimal lag-length choice in stable and unstable VAR models under situations of homoscedasticity and ARCH[END_REF].

The following restrictions on the deterministic parameters are considered

R l : µ 1 = ατ l Rl : µ 1 = 0 (2.2)
R s : µ 0 = ατ s and µ 1 = 0 Rs : µ 0 = µ 1 = 0.

Note that we do not have necessarily τ l = 0 for R l and τ s = 0 for R s or µ 0 = 0 for R l , Rl and R s , so that we have the relation Rs ⊂ R s ⊂ Rl ⊂ R l . Nested LR tests for the restrictions Rl , R s and Rs are proposed in [START_REF] Johansen | Likelihood-Based Inference in Cointegrated Vector Autoregressive Models[END_REF]. Now we discuss the consequences of the restrictions in (2.2) on the behaviour of (x t ). Suppose that α ′ ⊥ Γβ ⊥ has full rank, where

Γ = I d - p0-1 i=1 Γ i .
We also assume that the autoregressive polynomial

A(z) = (1 -z)I d -αβ ′ z -p0-1 i=1 Γ i (1 -z)z i is such that det A(z) = 0 implies that | z |> 1 or z = 1.
Under these assumptions, it follows from Granger's representation theorem [START_REF] Johansen | Likelihood-Based Inference in Cointegrated Vector Autoregressive Models[END_REF], p 49) that

x t = C t i=1 (ǫ i + µ 0 + µ 1 t) + D(L)(ǫ t + µ 0 + µ 1 t) + A = C t i=1 ǫ i + ρ 1 t + ρ 0 + Y t + A, (2.3) where C = β ⊥ (α ′ ⊥ Γβ ⊥ ) -1 α ′ ⊥
and L is the usual lag operator. The vector A depends on initial values and verifies β ′ A = 0. The process (Y t ) is linear and such that

y t = ∞ i=0 φ i ǫ t-i ,
where the power series

D(z) = ∞ i=0 φ i z i is convergent for | z |≤ 1 + κ for some κ > 0.
The vectors ρ 1 , ρ 0 are functions of the parameters in (2.1). If we suppose that R l hold, we may have ρ 1 = 0 and ρ 0 = 0 and in this case (β ′ x t ) can be composed of a stationary process plus a linear trend. We say in this case that (β ′ x t ) is trend stationary. If Rl hold, again we may have ρ 1 = 0 and ρ 0 = 0, but (β ′ x t ) is stationary. If R s hold we obtain ρ 1 = 0 but we still may have ρ 0 = 0, and in this case it is also allowed to have E(β ′ x t ) = 0. Finally the restriction Rs does not allow for any deterministic component for (β ′ x t ) and (x t ). The number of independent linear combinations β ′ x t which are such that the random walk behaviour is vanished is the cointegrating rank. In particular when r 0 = 0 there is no cointegration between the variables and model (2.1) is a vector autoregressive model for the process (∆x t ). From (2.3) the process (∆x t ) is stationary for all the restrictions, so that (x t ) is I(1). Note that we do not study the unrestricted case which produce nonstationary processes with quadratic trend since it is rarely faced in applied works. Then we see that the deterministic term and the cointegrating rank are important for the data analysis and forecasting purposes. However we will see in the next section that the specification of µ 0 , µ 1 and r 0 strongly depend on the choice of the autoregressive order.

Various tests for the cointegrating rank are available in the literature (see [START_REF] Hubrich | A review of systems cointegration tests[END_REF] for a review of such tests). The most commonly test for the cointegration rank is the LR test introduced by [START_REF] Johansen | Statistical analysis of cointegration vectors[END_REF][START_REF] Johansen | Estimation and hypothesis testing of cointegration vectors in gaussian vector autoregressive models[END_REF], which is shown to be asymptotically valid under quite general assumptions (see [START_REF] Rahbek | ARCH innovations and their impact on cointegration rank testing[END_REF] or [START_REF] Raïssi | Testing the cointegrating rank with uncorrelated but dependent errors[END_REF]). Using the same arguments of [START_REF] Raïssi | Testing the cointegrating rank with uncorrelated but dependent errors[END_REF] and following the proofs of corollary 11.2 and theorem 11.3 of [START_REF] Johansen | Likelihood-Based Inference in Cointegrated Vector Autoregressive Models[END_REF] it can be shown that the LR tests for the deterministic term are also asymptotically valid when the errors are dependent but uncorrelated.

Determining the autoregressive order

We first briefly describe the estimation procedure of the model (2.1). The reader is referred to [START_REF] Johansen | Likelihood-Based Inference in Cointegrated Vector Autoregressive Models[END_REF] for more details on the maximum likelihood estimation of VECM. If the errors are not assumed Gaussian, the quasi maximum likelihood is used. Let us assume that the observations x 1 , . . . , x n are available. Note that α and β are not identified in (2.1), and in case of R l and R s the parameters τ l and τ s are also not identified. In the sequel we suppose that these parameters are normalized in some appropriate way. Let us rewrite (2.1) as follow

∆x t = θξ t-1 + ǫ t , (3.1)
where

ξ t-1 = ((B ′ z t-1 ) ′ , ∆X ′ t-1 , . . . , ∆X ′ t-p0+1 ) ′ . We set B = (β ′ , τ l ) ′ and z t-1 = (x ′ t-1 , t) ′ if we suppose that R l hold but not Rl in (3.1). If Rl is hold but not R s we set B = β and z t-1 = x t-1 . We consider B = (β ′ , τ s ) ′ and z t-1 = (x ′ t-1 , 1) ′ if R s is hold but not Rs . Finally we set B = β and z t-1 = x t-1 if the deterministic
parameters are assumed to be equal to zero. The parameters in B can be estimated using Reduced Rank (RR) regression. The super-consistency of the parameters in B can be obtained under quite general assumptions on the error terms (see e.g. [START_REF] Rahbek | ARCH innovations and their impact on cointegration rank testing[END_REF] or [START_REF] Raïssi | Testing the cointegrating rank with uncorrelated but dependent errors[END_REF]). It is important to note that in the RR estimation method the process (ξ t ) is used so that the resulting estimators strongly depend on the fitted autoregressive order p > 0 which can be such that p = p 0 . As a consequence the computation of the LR tests statistics for the cointegrating rank and the deterministic terms depend on p. Therefore the choice of the autoregressive order is crucial for the cointegration analysis and the specification of the deterministic part of the model. Let

us set θ = [α, Γ 1 , . . . , Γ p0-1 , µ 0 ] if R l hold but not R s and θ = [α, Γ 1 , . . . , Γ p0-1 ] if R s hold.
Once the parameter B is replaced by its estimator in (3.1) we can compute a least squares estimator for θ. However it is clear that the estimator of the short run parameters we obtain depend on the specification of cointegrating rank and of the deterministic parameters. If we suppose that the VECM is well specified, it can be shown that these estimators are asymptotically normally distributed (see Brüggemann, [START_REF] Brüggeman | Residual autocorrelation testing for error correction models[END_REF] in the iid Gaussian case and Raïssi (2010) when the errors are dependent but uncorrelated). Now we introduce the tests used in the Box-Pierce methodology. Let us denote by ǫt the residuals obtained from the estimation stage. We define the residual autocovariances

Γǫ (h) := n -1 n t=h+1 ǫt ǫ′ t-h .
The commonly used portmanteau statistic based on the first m autocovariances is given by

Q = n m h=1 Tr Γǫ (h) ′ Γǫ (0) -1 Γǫ (h) Γǫ (0) -1 . (3.2)
This statistic corresponds to the generalization of the [START_REF] Box | Distribution of residual autocorrelations in autoregressive-integrated moving average time series models[END_REF] portmanteau statistic in the multivariate case proposed by [START_REF] Chitturi | Distribution of residual autocorrelations in multiple autoregressive schemes[END_REF]. One can alternatively use the [START_REF] Ljung | On measure of lack of fit in time series models[END_REF] portmanteau statistic proposed in the multivariate framework by [START_REF] Hosking | The multivariate portmanteau statistic[END_REF] for potential improvements when the errors are Gaussian. Nevertheless the use of the Ljung-Box statistic lead to the same general conclusion for our study, so that we focus on the Box-Pierce (BP hereafter) portmanteau test in the sequel. The tested hypotheses for checking the adequacy of the autoregressive order are given by

H 0 : E(ǫ t ǫ t-i ) = 0 vs H 1 : ∃ i such that E(ǫ t ǫ t-i ) = 0, i ∈ {. . . ,
(d 2 (m -p 0 + 1) -dr 0 ) distribution when m → ∞ as n → ∞.
Using this result they proposed a portmanteau test in the case of cointegrated variables.

This test will be denoted as BP S in the sequel. If we suppose that the errors are only uncorrelated and ǫ t 4+ν < ∞ for some ν > 0 and under other technical assumptions, it is shown in the framework of Rs in [START_REF] Raïssi | Autocorrelation based tests for vector error correction models with uncorrelated but nonindependent errors[END_REF], that the asymptotic distribution of the test statistic Q is a weighted sum of chi-squares. Starting from this result a portmanteau test which is robust to the possible presence of nonlinearities in the error terms is proposed. The extension to the other restrictions considered in this paper is straightforward and is described in an extended version of Raïssi (2010) available under request. This test will be denoted as BP W . In practice the RR estimators are derived for given 0 ≤ r < d which can be such that r = r 0 . It is clear from the estimation stage that the statistic Q depend on the specified deterministic term and the fitted cointegrating rank. The effects of such possible misspecifications will be studied in the next section.

To determine the number of short run dynamics to include in the model information criteria of the form

IC = log | Γǫ (0) | +(p + 1)c n /n.
are often used for a given specification of the deterministic parameters and cointegrating rank r. [START_REF] Paulsen | Order determination of multivariate autoregressive time series with unit roots[END_REF] showed that the information criterion provide a consistent estimator of p 0 if and only if c n is such that c n → ∞ and c n /n → 0. We have c n = d 2 log n for the BIC and c n = 2d 2 for the AIC, so that the BIC gives a consistent estimator of the lag length on the contrary to the AIC. Then we focus on the BIC in the sequel, noting that any other convergent information criterion will lead to the same general result. The BIC is computed for a set of possible values of p and the lag length which minimizes the BIC is selected. In general the residual empirical variance Γǫ (0) is often used for the computation of information criteria. Then similarly to the portmanteau test statistics, we see that well specifying the cointegrating rank and/or the deterministic term is important. We also note from the BIC expressions that this approach for choosing p is not intended to detect a serial correlation of the residuals or to control any error of first kind. However we can remark that the existence of the fourth moments is not needed to compute these information criteria. In the next section we compare the two approaches for the choice of p.

Monte Carlo experiments

The performances of the BIC and of the tests based on the residual autocorrelations for selecting the autoregressive order is studied. Using these results we propose strategies which combine the use of the BP methodology and information criteria for selecting the lag length, taking into account for possibly misspecified deterministic terms or fitted cointegrating rank. In a first time the framework of the experiments is described. In each of the experiments we simulated N = 1000 independent trajectories of lengths n = 100 and n = 1000 using a Data Generating Process (DGP) inspired by the canonical form of Toda (1994)

x t = ψI r0 0 0 I d-r0 x t-1 + Γ 01 ∆x t-1 + µ 0 + µ 1 t + ǫ t ,
where ψ is a scalar and d = 4. The parameters are taken ψ = 0.5, Γ 01 = -0.3I 4 and the cointegrating rank is r 0 = 2. We set µ 1 = (-0.05, -0.05, 0, 0) ′ and µ 0 = (0.1, 0.1, 0.1, 0.1) ′ if R l is hold but not Rl . The results for the case of simulated processes such that R s is hold but not Rs are not displayed since they are similar to the case µ 0 = µ 1 = 0. In the sake of conciseness the results for the restriction Rl is hold but not R s are not displayed since they lead to similar conclusions to the cases presented here. This kind of DGP is used in numerous studies (e.g. [START_REF] Hubrich | A review of systems cointegration tests[END_REF] or [START_REF] Demetrescu | Testing for the cointegrating rank of a vector autoregressive process with uncertain deterministic trend term[END_REF]). In addition the results concerning the information criteria are corroborated by previous studies. In the experiments the cointegrating rank and the deterministic terms are not necessarily well fitted, so that the consequences of such situations is analyzed. Since the number of parameters is high, a relatively large samples are taken in our experiments. Two cases are considered for the error process in the sequel. In case of iid innovations we take ǫ t ∼ N (0, I d ). To illustrate the case where the errors are dependent, we use the following ARCH model with constant correlation of Jeantheau (1998)

   ǫ 1t . . . ǫ dt    =    σ 1t 0 0 0 . . . 0 0 0 σ dt       η 1t . . . η dt    (4.1)
where

   σ 2 1t . . . σ 2 dt    =    0.1 . . . 0.1    +     0.2 0 0 0.1 0.1 0.2 0 0.1 0.1 0 0.2 0.1 0.1 0 0 0.2        ǫ 2 1t-1 . . . ǫ 2 dt-1    ,
and the process η t = (η 1t , . . . , η dt ) ′ is iid, such that η t ∼ N (0, I d ).

The autoregressive order selection strategy correspond to the case where the practitioner have no prior knowledge on p 0 . More precisely for the portmanteau tests p = 1, 2, . . . are successively tested until the null hypothesis is not rejected. However we assume that when the autoregressive order adequacy is rejected for all p ≤ 5, the model specification is suspected to be not reliable and one stop the procedure. Indeed when p takes large values the number of parameters become high beside the available observations. In such case the analysis of the linear dynamics is not reliable. In our experiments we used the asymptotic nominal level 5% for the different portmanteau tests considered in this paper. In the sequel note that if the portmanteau tests perform well, the selected p = p 0 is often around 95% since this selection procedure is based on the control of the error of first kind. Note that the BP W test is available for fixed m, while for the BP S it is assumed that m → ∞ as T → ∞. Therefore m = 5 is taken for the BP W test for T = 100 and T = 1000, while we take m = 5 when T = 100 and m = 15 when T = 1000 for the BP S test. The autoregressive order which minimize the BIC in the set {1, . . . , 5} is selected. The results are given in Tables 123456789101112.

We first analyze the results for small samples (T = 100) given in Tables 123456. It appears that the performances of the BP W test is rather disappointing. This can be explained by the fact that about 40 parameters are estimated when p = 2 and thus the number of observations is relatively small for the elaborated BP W test. In such situation the more simple BP S test provide satisfactory results. It seems that the results for the portmanteau tests are not much affected when the deterministic terms or the cointegrating are misspecified in small samples. We only remark that when a deterministic trend present in the data is not taken into account by the practitioner (Tables 456), the BP S test tends to select a too large lag length. From Tables 1-6 we see that the BIC selects a too small p. When r < r 0 is taken the BIC is likely to choose a larger autoregressive order. Therefore it seems that the use of the parsimonious BIC when the sample is small can be misleading for the analysis of VECM. Now we consider the results for large samples (T = 1000). It is found from Tables 7-12 that the BP tests perform well when the long run relations and the deterministic terms are well specified. However we note that when r < r 0 is taken, the BP tests tend to select too large autoregressive orders. If r < r 0 the estimators of the short run parameters are biased and the residual autocorrelations may appear too large as the sample increase even if p = p 0 . When r > r 0 is taken we note that the BP tests are not much affected. A possible explanation is that in this case the effect of non stationary components introduced in the model are limited by the fact that the corresponding adjustment parameters are close to zero. It is important to note that overspecified deterministic terms seems not entail a significant loss of efficiency in selecting the autoregressive order using the BP methodology. However from Tables 10-12 we see that the BP tests are clearly likely to select a too large autoregressive order when a deterministic trend is present in the observed process and is not taken into account in the model. Similarly we also found from Tables 10-12 that when deterministic trends are spuriously assumed to not enter the cointegrating relations the BP tests select too large p. When the BP tests are used and in case of doubt, it seems preferable to consider the restriction R l only for fitting the autoregressive order. In such situations the unnecessary parameters introduced in the model are few. Actually we noted no major loss of efficiency and considering such restriction allow to avoid misspecified autoregressive order. Of course as pointed out by several studies (e.g. [START_REF] Johansen | The role of the constant and linear terms in cointegration analysis of non-stationary variables[END_REF]) it is important to remove any misspecification of the deterministic term when determining the cointegrating rank. Finally it emerges that in general the BP W test is not affected by the presence of dependent errors. In some cases the BP S test tends to select too large p when the errors follow an ARCH model. This can be explained by the fact that the BP S is not intended to take into account such situations. In accordance with the theoretical, we found that the BIC choose successfully the optimal lag length when the sample is large in both standard and dependent cases. It also appears that the selected autoregressive orders do not depend much on the specification of the deterministic terms or on the fitted cointegrating rank. Nevertheless it is important to note that the use of information criteria for fitting the lag length can be quite misleading even when the sample is large. For instance we give the autoregressive order selected by the AIC criterion for the case given in Table 10. Despite the VECM is well specified, the AIC select a too large p. Results not reported here show that the AIC almost always choose a lag length p = 5 in many of the situations presented here. This is not surprising since for the AIC the more complicated models are not so penalized. Thereby this confirms that the BIC perform better than the AIC for the lag selection (see e.g. [START_REF] Cheung | Finite-sample sizes of Johansen's likelihood ratio tests for cointegration[END_REF] or Ho and Sorensen (1996)).

Conclusion

In this paper we studied the problem of the fit of the autoregressive order by comparing tools based on the Box-Pierce methodology recently proposed in the literature and information criteria commonly used in practice. Fitting an adequate autoregressive order is important for the analysis of VECM, but this task is carried out using information criteria only in general. It emerges from our study that underspecified deterministic trend or cointegrating rank can be quite misleading when using the Box-Pierce methodology for the choice of the autoregressive order p. We also found that the selection of p using information criteria does not depend much on the specification of the other parts of the model. However the use of these information criteria can be quite misleading in some cases. In general we recommend to use the information criteria together with the portmanteau tests when choosing the lag length of a VECM. When the selected autoregressive orders are close and small, then one should take the p proposed by the portmanteau tests. It is advisable to choose the p selected by the portmanteau tests since they are able to detect residual autocorrelations which can strongly affect the cointegration analysis. If the lag lengths selected by the two approaches presented in this paper are both large, it is important that the chosen information criterion is parsimonious to conclude on the reliability of the selected lag length. If the p retained by an information criterion is much larger than the p retained by the portmanteau tests, then one can suspect that the more complicated model is not penalized enough by the information criterion. On the other hand if the p selected using the portmanteau tests is much larger than the p selected by using an information criterion, one can suspect that the deterministic trend is underspecified or the fitted cointegrating rank is smaller the the true cointegrating rank.

Tables Table 1: Frequency (in %) of selected lag length using the portmanteau tests. The simulated processes are of length T = 100, such that Rs hold and r0 = 2. BP W 47.1 51.8 0.9 0.1 0.1 53.7 45.6 0.6 0.0 0.1 BP S 2.4 90.9 6.7 0.0 0.0 2.2 89.2 8.6 0.0 0.0 R l but not Rl BIC 98.8 1.2 0.0 0.0 0.0 98.0 2.0 0.0 0.0 0.0 Rl but not R s BIC 97.9 2.1 0.0 0.0 0.0 96.2 3.8 0.0 0.0 0.0 R s but not Rs BIC 98.1 1.9 0.0 0.0 0.0 96.5 3.5 0.0 0.0 0.0 µ 0 = µ 1 = 0 BIC 97.1 2.9 0.0 0.0 0.0 95.8 4.2 0.0 0.0 0.0

Case

RR

  l but not Rl BIC 96.6 3.4 0.0 0.0 0.0 95.1 4.9 0.0 0.0 0.0 Rl but not R s BIC 96.4 3.6 0.0 0.0 0.0 94.9 5.1 0.0 0.0 0.0 R s but not Rs BIC 97.0 3.0 0.0 0.0 0.0 95.3 4.7 0.0 0.0 0.0 µ 0 = µ 1 = 0 BIC 96.4 3.6 0.0 0.0 0.0 94.9 5.1 0.0 0.0 0.0Table 2: Frequency (in %) of selected lag length using the portmanteau tests. The simulated processes are of length T = 100, such that Rs hold and r0 = 2. l but not Rl BIC 57.7 43.2 0.0 0.0 0.0 56.6 43.4 0.0 0.0 0.0 Rl but not R s BIC 56.8 43.2 0.0 0.0 0.0 55.2 44.8 0.0 0.0 0.0 R s but not Rs BIC 59.4Table3: Frequency (in %) of selected lag length using the portmanteau tests. The simulated processes are of length T = 100, such that Rs hold and r0 = 2.

RR

  l but not Rl BIC 48.6 51.4 0.0 0.0 0.0 44.2 55.8 0.0 0.0 0.0 Rl but not R s BIC 45.1 54.9 0.0 0.0 0.0 43.8 56.2 0.0 0.0 0.0 R s but not Rs BIC 55.5 44.5 0.0 0.0 0.0 57.1 42.9 0.0 0.0 0.0 µ 0 = µ 1 = 0 BIC 54.3 45.7 0.0 0.0 0.0 56.9 43.1 0.0 0.0 0.0Table6: Frequency (in %) of selected lag length using the portmanteau tests. The simulated processes are of length T = 100, such that R l hold but not Rl and r0 = 2. Rl but not R s BIC 71.7 28.3 0.0 0.0 0.0 66.0 34.0 0.0 0.0 0.0 R s but not Rs BIC 73.4 26.6Table7: Frequency (in %) of selected lag length using the portmanteau tests. The simulated processes are of length T = 1000, such that Rs hold and r0 = 2. l but not Rl BIC 0.0 100.0 0.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 Rl but not R s BIC 0.0 100.0 0.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 R s but not Rs BIC 0.0 100.0 0.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 µ 0 = µ 1 = 0 BIC 0.0 100.0 0.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0

R l but not Rl BIC 0

 0

  -1, 0, 1, . . . } are commonly considered for the portmanteau tests, although the test statistics are only based on the m first residual autocorrelations. Assuming that the error process is iid Gaussian and such that ǫ t 4 < ∞, Brüggemann,[START_REF] Brüggeman | Residual autocorrelation testing for error correction models[END_REF] showed that the asymptotic distribution of the statistic in (3.2) can be approximated by a χ 2

Table 4 :

 4 Frequency (in %) of selected lag length using the portmanteau tests. The simulated processes are of length T = 100, such that R l hold but not Rl and r0 = 2.

		Case			iid errors				ARCH errors
	r = 2	p	1	2	3	4	5	1	2	3	4	5
	R l but not Rl Rl but not R s	BP W l BP S l BP W l BP S l	0.2 98.5 1.1 0.1 0.1 0.0 94.4 5.6 0.0 0.0 0.2 97.4 1.7 0.2 0.5 0.0 90.9 9.1 0.0 0.0	2.7 94.3 2.1 0.7 0.2 0.0 89.4 10.5 0.1 0.0 2.7 94.6 1.7 0.3 0.7 0.0 87.6 12.3 0.1 0.0
	R s but not Rs	BP W 21.3 76.7 1.6 0.1 0.3 26.3 71.6 1.5 0.3 0.3 BP S 0.5 92.6 6.9 0.0 0.0 0.3 82.6 16.9 0.2 0.0
	µ 0 = µ 1 = 0	BP W 16.6 81.6 1.4 0.2 0.2 21.7 76.4 1.4 0.3 0.2
		BP S	0.4 91.6 7.9 0.1 0.0	0.2 83.5 16.1 0.2 0.0
	R l but not Rl	BIC 89.

Table 5 :

 5 Frequency (in %) of selected lag length using the portmanteau tests. The simulated processes are of length T = 100, such that R l hold but not Rl and r0 = 2.

		Case			iid errors				ARCH errors
	r = 1	p	1	2	3	4	5	1	2	3	4	5
	R l but not Rl Rl but	BP W BP S BP W	0.2 97.4 1.8 0.3 0.3 0.0 94.1 5.9 0.0 0.0 0.2 96.8 2.1 0.4 0.5	2.7 95.1 1.8 0.4 0.0 0.0 91.8 8.2 0.0 0.0 2.7 95.0 1.6 0.3 0.4
	not R s

Table 8 :

 8 Frequency (in %) of selected lag length using the portmanteau tests. The simulated processes are of length T = 1000, such that Rs hold and r0 = 2.

		Case			iid errors				ARCH errors
	r = 1	p	1	2	3	4	5	1	2	3	4	5
	R l but not Rl Rl but	1.7 0.0 35.8 39.4 18.5 6.3 2.5 3.3 92.5 0.0 BP W 0.0 0.0 32.6 41.7 19.1 6.6 3.7 2.2 3.0 91.1 BP S BP W 0.0 1.9 2.4 3.3 92.4 0.0 3.6 2.3 3.0 91.1
	not R s	BP S	0.0 35.4 39.7 18.5 6.4	0.0 32.5 42.0 19.0 6.5
	R s but not Rs	1.9 0.0 36.0 39.6 17.9 6.5 2.4 3.3 92.4 0.0 BP W 0.0 BP S 0.0 32.7 42.3 18.8 6.2 3.7 2.2 3.0 91.1
	µ 0 = µ 1 = 0	BP W 0.0	1.6	2.6	3.1 92.7 0.0	3.6	2.1	3.2 91.1
		BP S	0.0 35.8 39.7 18.2 6.3	0.0 32.8 42.5 18.2 6.5

Table 9: Frequency (in %) of selected lag length using the portmanteau tests. The simulated processes are of length T = 1000, such that Rs hold and r0 = 2.

Case

iid errors ARCH errors r = 3 p 1 2 3 4 5 1 2 3 4 5 R l but BP W 0.0 95.7 1.4 0.9 2.0 0.0 96.5 1.4 0.7 1.4 not Rl BP S 0.0 94.3 4.6 1.1 0.0 0.0 93.2 5.7 1.0 0.1 Rl but BP W 0.0 95.8 1.4 0.7 2.1 0.0 96.2 1.2 1.0 1.6 not R s BP S 0.0 94.4 4.5 1.0 0.1 0.0 92.7 6.2 1.0 0.1 R s but BP W 0.0 95.8 1.4 0.7 2.1 0.0 96.0 1.5 0.9 1.6 not Rs BP S 0.0 94.4 4.5 1.1 0.0 0.0 92.7 6.2 1.0 0.1 µ 0 = µ 1 = 0 BP W 0.0 95.8 1.4 0.8 2.0 0.0 96.1 1.4 0.8 1.7 BP S 0.0 94.0 5.0 1.0 0.0 0.0 93.0 5.9 1.0 0.1 R l but not Rl BIC 0.0 100.0 0.0 0.0 0.0 0.1 99.9 0.0 0.0 0.0 Rl but not R s BIC 0.0 100.0 0.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 R s but not Rs BIC 0.0 100.0 0.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 µ 0 = µ 1 = 0 BIC 0.0 100.0 0.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0