Thomas Barré
email: thomas_barre@yahoo.fr

Price expectations and price dynamics: the case of the rice sector in developing Asia

Keywords: JEL Classification: C23, D41, D84, E31, O13, Q11

Uncertainty is a crucial issue for producers who must make input decisions without knowing prices and without perfect knowledge of realized output. In this context, price expectations strongly determine the production choices and market prices that result from market-clearing conditions. This study analyzed the role that price expectations play in price dynamics, developing a theoretical model of trade in varieties following Armington (1969) and augmented with yield and price uncertainty to highlight several main determinants of domestic producer prices, including exchange rates, proximity to world markets, input prices, natural disasters, and producers' expectations. An econometric estimation of the rice sector, using a panel of 13 developing Asian countries during 1965-2003, confirmed that expectations count, with a 1% increase in the expected price resulting in a 1.18% decrease in the market price. A simulation exercise based on these empirical results demonstrated that forecasting errors are large. Specifically, Asian rice farmers have a 50% chance of making prediction errors of 10% or more on the final market price. This high error rate suggests the need for developing ways of sharing information, such as radio programs dedicated to agricultural producers or the introduction of futures markets, to stabilize agricultural incomes.

The food crisis of 2007-2008 and the recent surge in food prices have renewed the interest of economists in employing policy tools or instruments appropriate for managing agricultural price volatility. As proposed by [START_REF] Galtier | How to manage food price instability in developing countries? Working Paper MOISA[END_REF], these policy tools can be classified into four categories based on their objective (price stabilization or reduction in the effect of price instability) and their -mode of governance‖ (public or market based). The problem often pointed with price stabilization is that it can distort the informational content of prices. If the state intervenes to stabilize priceseither through public buffer stocks or trade restrictionsproducers are incentivized to either overproduce when the state guarantees a floor price to producers or under produce when it sets a ceiling price for consumers. This policy, however, disrupts the negative correlation between prices and production levels and can cripple the government revenues, as was the case in [2007][2008]. On the other hand, market-based instruments like futures hedging coupled with transfers to the most vulnerable populations unable to use market-based instruments appears to be the best way to protect the poor against price volatility without affecting prices. But the lack of acceptance of such instruments aimed at protecting peasants against price volatility by developing countries [START_REF] Crmg | The International Task Force on Commodity Risk Management in Developing Countries: Activities, Findings and the Way Forward[END_REF] has led researchers to focus on market based instruments that can stabilize prices without affecting their informational content. These policies are supposed to contribute towards a good match between supply and demand and stabilize prices thanks to improvements in spatial and temporal integration of markets by the development of transport networks, storage facilities and quality standards. Finally, the recent surge in food prices has demonstrated that while state intervention remains necessary to avoid domestic food crises and protect the most vulnerable populations, it can also induce hoarding behavior among consumers and speculators, putting additional pressure on prices [START_REF] Timmer | Did Speculation Affect World Rice Prices? Agricultural and Development Economics Division of the Food and Agriculture Organization of the United Nations (FAO -ESA)[END_REF], providing a new argument against state intervention.

As the usefulness of price management policies depends on the sources of price dynamics, the core issue in the debate on the use of the appropriate instruments to manage price volatility is to determine the forces behind the formation of prices. [START_REF] Galtier | How to manage food price instability in developing countries? Working Paper MOISA[END_REF] distinguishes three main forces that determine price changes: The first cause of agricultural price volatility is -natural instability‖, which comes from natural hazards that affect production. When a natural disaster occurs, production decreases, resulting in excess demand with respect to supply and an increase in market prices. A second source of instability of domestic prices is imported, that is, it is caused by changes in parity prices and substitutions in consumption baskets by world customers. Trade policies and exchange rates are key factors in this form of instability. A third source of instability is the endogenous fluctuations generated by the instability of market players' price expectations.

In this case, -domestic prices can be unstable without any movement in market fundamentals (domestic supply, demand curves, and terms of trade)‖. For example, farmers increase their production level if price expectations increase, which in turn generates a drop in prices because supply exceeds demand. This article focuses on the role of price expectations and develops a model for agricultural commodities price formation based on an economic geography framework (following Armington 1969) augmented with price and yield uncertainty. The model was tested on an unbalanced panel dataset of producer prices for rice in 13 Asian countries. The theoretical model appears relevant as demonstrated by the estimation results. First, estimated structural parameters corresponded to those found in the literature and to prior knowledge, with an Armington elasticity close to 6 and returns to scale approximately equal to 0.9. Second, a simulation exercise determined that the quality of price expectations was so that farmers have a 50% chance of making a forecasting error being greater than 10% of the market price and showed that the model can reproduce some of the main characteristics of price changes (volatility, skewness, and kurtosis). The analysis of the data also showed that price expectations were a crucial determinant of price formation. Given the extent of forecasting errors, improvements in the informational network should be at the top of the reform agenda. In particular, the development of futures markets in the rice sector, combined with investments in radio programs dedicated to farmers, could address the poor price expectations issue and help in price stabilization, which is an important condition for sustainable growth [START_REF] Ramey | Cross-Country Evidence on the Link Between Volatility and Growth[END_REF].

The paper is organized as follows: Section 1 briefly presents the literature on agricultural commodity price formation; Section 2 describes the theoretical model and estimation method; Section 3 presents the data and empirical estimation results; and Section 4 concludes.

Taxonomy of price dynamics explanations

As noted by [START_REF] Gouel | Agricultural price instability: a survey of competing explanations and remedies[END_REF], most of the literature on agricultural price instability is characterized by a dichotomy between endogenous and exogenous explanations for price dynamics. In general, exogenous explanations identify harvest shocks as the cause of price fluctuations, and endogenous explanations identify errors in expectations as the cause. Other researchers have focused on monetary explanations for price changes [START_REF] Frankel | Expectations and Commodity Price Dynamics: The Overshooting Model[END_REF] or the role of exchange rates and the terms of trade [START_REF] Liefert | The Transmission of Exchange Rate Changes to Agricultural Prices[END_REF]).

The original model of endogenous fluctuations can be traced back to [START_REF] Ezekiel | The Cobweb Theorem[END_REF] and his cobweb theorem, in which it was assumed that farmers formed their price expectation on the basis of the observed actual price, . According to this theorem, when the actual price is high, farmers anticipate future prices to be high and increase their production level. This increase in the quantity harvested decreases the price. In the next period, given the low actual price, farmers anticipate a low future price and so will plan to produce less, which will in turn increase the price of the harvested good, and so on. The model was later extended by [START_REF] Nerlove | Adaptive Expectations and Cobweb Phenomena[END_REF], who proposed adaptive expectations, meaning that producers build their expectations as a function of the last expected price and last period forecasting error. One important implication of these cobweb models, even in their modern forms that include risk aversion [START_REF] Boussard | When risk generates chaos[END_REF] or non-linear curves [START_REF] Hommes | Dynamics of the cobweb model with adaptive expectations and nonlinear supply and demand[END_REF], is that they always predict negative first order autocorrelation of prices, while empirical evidence strongly supports positive autocorrelation (Deaton and Laroque 1996).

Another important drawback of backwards looking expectations is that, in the case of exponentially increasing inflation, producers systematically underestimate the future price. These systematic errors imply that backwards-looking expectations waste information and generate irrational behaviors. This is one of the main reasons for the development of the rational expectation model [START_REF] Muth | Rational Expectations and the Theory of Price Movements[END_REF], in which producers base their expectations on all the available information at the time of forecasting so that they cannot make systematic prediction errors:

, where is an i.i.d prediction error term. However, the introduction of rational expectations in a simple linear model reduces price dynamics to simple fluctuations around a steady state. Therefore, the simple rational expectations model is unable to reproduce the high degree of autocorrelation observed in commodity prices. This limitation led to the development of the workhorse in the literature on commodity price dynamics: the competitive storage model.

As noted by [START_REF] Muth | Rational Expectations and the Theory of Price Movements[END_REF], the introduction of storage in the linear rational expectation model generates positive autocorrelation in prices: -Speculation smoothes shocks over several periods, so the effect of one shock is spread across subsequent periods, causing positive autocorrelation‖ [START_REF] Gouel | Agricultural price instability: a survey of competing explanations and remedies[END_REF]. [START_REF] Deaton | On the Behaviour of Commodity Prices[END_REF] demonstrated that the introduction of competitive storage also affects the distribution of prices, increasing its kurtosis and allowing for positive skewness, two of the major characteristics of commodity prices. Despite the difficulties in estimating with this model due to its non-linear components, including structural breaks and the lack of data on inventories, many developments around the competitive storage model have been proposed in the past 20 years [START_REF] Gouel | Agricultural price instability: a survey of competing explanations and remedies[END_REF]. A common feature of the models based on competitive storage and rational expectations is that natural disasters are the only source of uncertainty and the only cause of errors in price expectations, assuming no endogenous explanations of price dynamics. [START_REF] Boussard | Storage and the Volatility of Agricultural Prices: a Model of Endogenous Fluctuations[END_REF] attempted to fill this gap by introducing adaptive expectations in a competitive storage model. However, their simulation results are not more accurate than the standard framework of [START_REF] Deaton | On the Behaviour of Commodity Prices[END_REF] since price autocorrelation in simulated series remains rather low when compared to actual data.

In all of these studies, natural disasters are the only cause for uncertainty and speculation in inventories. [START_REF] Frankel | Expectations and Commodity Price Dynamics: The Overshooting Model[END_REF] proposed an additional mechanism, the overshooting hypothesis, based on monetary shocks, that is summarized by the author as:

"A monetary contraction temporarily raises the real interest rate (whether via a rise in the nominal interest rate, a fall in expected inflation, or both). Real commodity prices fall. How far?

Until commodities are widely considered "undervalued" --so undervalued that there is an expectation of future appreciation (together with other advantages of holding inventories, namely the "convenience yield") that is sufficient to offset the higher interest rate (and other costs of carrying inventories: storage costs plus any risk premium). Only then are firms willing to hold the inventories despite the high carrying cost. In the long run, the general price level adjusts to the change in the money supply. As a result, the real money supply, real interest rate, and real commodity price eventually return to where they were."

In summary, the overshooting hypothesis feeds the competitive storage model with, apart from yield uncertainty, an additional source for speculation from unanticipated monetary shocks, which can cause price spikes and falls.

Apart from these explanations based on speculative storage, Deaton and Laroque (2003) showed that it is also possible to represent positive autocorrelation, skewness, and kurtosis of observed data series with a rational expectations model only made of producers and consumers, i.e. assuming away competitive storage. Deaton and Laroque (2003) proposed a modified version of the [START_REF] Lewis | Economic Development with Unlimited Supplies of Labour[END_REF] model that assumes a linear stochastic demand function (increasing with income) and a stochastic supply (assuming non-maximizing behavior) that equals its previous value corrected for the excess of the current price over the marginal cost of production plus a possibly autocorrelated supply shock. Thus, the standard -production lag‖ approach, according to which producers plan production before prices are revealed, is assumed away. Using this setting, -price behavior comes from the action of an integrated (trending) demand process against a supply function that is infinitely elastic in the long run but not the short run.‖ However, empirical evidence is insufficient to provide -any direct statistical support for the model‖ (Deaton and Laroque, 2003).

According to [START_REF] Galtier | How to manage food price instability in developing countries? Working Paper MOISA[END_REF], endogenous and exogenous explanations of price dynamics are accompanied by -imported instability‖, which encompasses many situations from weather shocks abroad to trade restrictions that can be summarized into a terms-of-trade effect. First, a supply shock in a foreign country can be partially transmitted to domestic prices because consumers substitute goods when relative prices change. Second, when the local currency is depreciated with respect to the reference unit, imported goods become more expensive and exports become cheaper. Consequently, the demand for locally produced goods increases, which in turn puts pressure on domestic prices if supply is not perfectly elastic or lags in response to price changes such as frequently occur in the rice sector because of production delays. Finally, when trade restrictions are strengthened domestically, foreign demand for local goods and local demand for goods produced abroad decrease, changing domestic prices.

The present study follows Deaton and Laroque (2003) in assuming away speculative storage and concentrates on the role of price expectations that is emphasized by the literature on endogenous expectations, in a rational expectations framework. A model of international trade in commodities in the spirit of [START_REF] Armington | The Geographic Pattern of Trade and the Effects of Price Changes[END_REF] is developed and augmented with yield and price uncertainty, so that it allows for the 3 sources of price uncertainty listed by [START_REF] Galtier | How to manage food price instability in developing countries? Working Paper MOISA[END_REF], i.e.

expectations, production shocks, and trade.

Theoretical framework and estimation strategy

The theoretical framework is based on an international trade model in commodities that follows [START_REF] Armington | The Geographic Pattern of Trade and the Effects of Price Changes[END_REF]. Each country produces its own unique variety of each good. Within each country, a unique, representative consumer presents preferences over the whole set of goods available in the world and makes consumption choices once prices have been revealed.

Consumers all share the same preferences. Production of a country-specific variety of a specified good is carried out by a unique competitive firm (no market power). This firm faces a production lag problem in that it plans its production quantities before prices have been revealed. The firm also faces different sources of uncertainty, including yield uncertainty and price uncertainty, which is partly the result of yield uncertainty as well as unpredictable demand.

I start by presenting consumer behavior and the resulting demand function addressed by each producer once prices have been revealed. The maximization problem of producers is then revealed, and the equilibrium market price equation is detailed together with the estimation strategy.

Consumer behavior

The consumer in region maximizes his utility at time , , which is represented by a Constant Elasticity of Substitution function over every good available in the world, which is made up of regions, where the region produces goods. represents the quantity of the variety of good produced in region that is consumed by the representative consumer of region at time . represents the weight attributed to this variety in the consumer utility function. represents the constant elasticity of substitution. When tends to 1, goods are perfect complements, and when tends to infinity, goods are perfect substitutes.

(1)

Consumers maximize their utility with respect to a budget constraint, where is the available income of the consumer in region at time . represents the price of the variety of commodity produced in region in local currency. represents the exchange rate of region currency with the US dollar at time ( units of region s local currency equals one USD).

represents an iceberg transport cost from region to region at time . Solving the cost minimization program of this consumer, we obtain:

(2)

Which gives:

(3)

Consumer demand at time of the variety of commodity produced in region is finally:

(4)

This demand function decreases with the commodity price, transport cost, and exchange rate of country with the USD, and increases with the share parameter , income, price index , and region exchange rate.

Given iceberg transport costs, the demand addressed by world consumers to the producer of good in country is:

(5) and total exports from country to country at time is:

(6)

This gravity trade equation must be estimated for each year using bilateral trade flows, where is an importer-year fixed effect, is an exporter-year fixed effect, and is a bilateral iceberg transport cost. In logarithmic terms, the equation is:

(7)

This equation has been estimated by [START_REF] Head | Gravity, market potential and economic development[END_REF] for all countries in the world with available trade data during 1960-2003. [START_REF] Head | Gravity, market potential and economic development[END_REF] defined and computed -Market Potential‖ data (a measure of proximity to world markets) for each country in each year:

(8)

Our world demand function reduces to:

Where is estimated by [START_REF] Head | Gravity, market potential and economic development[END_REF], and prices and exchange rates data are available from standard sources, as discussed in the next section.

Producer Behavior

A competitive firm maximizes expected profits given a stochastic Cobb-Douglas production function, which depends on the use of a composite input bought in at the price and unpredictable production shocks . This composite input might include, among other production factors, seeds, labor, arable land, fertilizers, energy, machinery, pesticides, and herbicides. A return to scale parameter, , is attached to and is expected to be inferior but close (or equal) to one in order to be consistent with the non-increasing (or constant) return to scale hypothesis often cited in analyses of agricultural production functions [START_REF] Bardhan | Size, Productivity, and Returns to Scale: An Analysis of Farm-Level Data in Indian Agriculture[END_REF]Townsend, Kirsten, and Vink 1998). Ex post profits equal the revealed price multiplied by the amount of harvest sold minus the cost of inputs bought in , .

First order conditions of maximization with respect to , give:

(11)

After some substitutions, the equation is:

(12) This is the supply function of the producer of good in region . It depends on input prices, the amount of past harvest kept for seeding, price expectations, and harvest shock. We assume to be constant [i.e., we assume that is constant].

Market equilibrium

The market clearing condition is , which gives the equilibrium price equation, taking logs and first differentiating:

(13)

This equation highlights five main determinants of commodity prices. First, in line with the literature, exchange rates play a crucial role on exchange rate pass-through to commodity prices. When the local currency is depreciated increases) with respect to the reference (the US dollar), then world demand for national products increases and domestic prices become more expensive. However, the amount of this exchange rate pass-through is different from (but close to) one since the exchange rate also determines the market access variable, which measures the proximity of the country to -world income‖. As expected, an increase in world wealth is transmitted to domestic prices, and the elasticity of prices with respect to world income is the inverse of the constant elasticity of substitution characterizing the consumer's utility function .

Empirical evidence regarding the value of in the literature on international trade will help us in determining whether our estimations are consistent with our framework or not. Also, our framework predicts that a 1% increase in the country's Market Access is equivalent (in terms of price change) to a 1% change in the production shock . Input prices play an important role in this equation. If returns to scale, are close to (but less than) one, the elasticity of producer price with respect to input prices can tend towards infinity. Finally, equation ( 13) quantifies the role of price expectations in price dynamics. This negative correlation comes from the interaction between the behavior of producers, who increase their production level when they expect prices to rise, and market conditions, which following an increase in production induce a fall in prices due to excess supply with respect to demand. As in the case of input prices, price forecasts are expected to be an important determinant of price dynamics, since price elasticity with respect to expectations can tend to (minus) infinity.

This equation cannot be estimated in its actual form since we don't observe price expectations. Therefore, the core issue at this step is to define a price forecasting rule. I adopt rational expectations, which have two virtues in this setting. First, it is the most -optimistic‖ framework since it assumes that producers don't make systematic forecasting errors, as would be the case with naïve or adaptive expectations. Second, it is the only framework that doesn't require dynamic panel data methods, which are difficult to handle given that available data have neither large cross-section nor important time dimensions (more details in the following section).

Following [START_REF] Lovell | Tests of the Rational Expectations Hypothesis[END_REF], I define the rational expectations' forecast error, so that

, where is identically and independently distributed, and replace the expected price by its expression as a function of the true price and the expectation error in equation ( 13). I obtain an expression for the price equation that has been purged from the presence of unobserved price expectations:

(14) Equation ( 14), which can be estimated by Ordinary Least Squares, draws attention to the normative definition of price expectations. Price expectations are thought to be -good‖ if forecasting errors , are close to unity. One of the objectives of the present study is to measure -how good‖ price expectations are.

Estimating equation ( 14) by Ordinary Least Squares (OLS) would usually be considered misleading due to endogeneity biases. Indeed, forecasting errors, , are correlated with unanticipated changes in the exchange rate , market access , and production shocks However, equation ( 13) gives a structural measure of this endogeneity bias so that OLS is efficient when estimating equation ( 14). Coefficients of equation ( 13) can be retrieved using -biased‖ coefficients of equation ( 14), even if we cannot observe the expected price. Also, OLS coefficients attached to and permit computation of structural parameters , the consumer's elasticity of substitution, and , the extent of the returns to scale, which can be compared to prior knowledge to check for consistency of empirical results obtained with the theoretical framework.

Data and Empirical results

While most analyses of commodity prices focus on international prices, the framework proposed here allows for the study of country-specific producer data, which is better for determining the implications of policies on producers. Given that our theoretical model assumed the same production technology for all producers of the same good all over the world, we focused on a single commodity, rice, in a geographically and economically narrowed region, developing Asia. We collected producer price (IRRI website) and quantities harvested (FAOSTAT) for rice in 13 Asian countries between 1965 and 20031 (Table 1). This choice is motivated by the need for a homogenous region of study because our theoretical model assumes a common production technology for every country. Rice is an appropriate commodity for this kind of assumption given that there are few differences in the production systems used to cultivate rice in this region (except in Japan where production is highly mechanized and which is excluded from the sample). 1971-2003Nepal 1965-2003Cambodia 1993-2003Pakistan 1966-2003China 1966-2003Philippines 1965-2003India 1965-2003Sri Lanka 1965-2003Indonesia 1967-2003Thailand 1965-2003Laos 1990-2003Vietnam 1991-2003Malaysia 1973-2003 Exchange rates data were obtained from IMF-IFS. Market potentials calculated by [START_REF] Head | Gravity, market potential and economic development[END_REF] were obtained from the CEPII website. Equation ( 14) depended on a composite input price , which was not observed. We used available data on input prices (fertilizers, energy, and wages) as a proxy for this variable. However, input price data were scarce, especially for fertilizer and wages. To address this issue, we used agricultural gross domestic product per capita, computed as the agricultural GDP in local currency divided by the agricultural labor force (obtained from IRRI), as a proxy and filled in gaps in fertilizer prices (urea prices were obtained from IRRI) using data from other countries. We converted fertilizer prices into US dollars, computed a year specific mean, and filled gaps by converting those means into local currencies. Energy prices were approximated by the petroleum spot (US$/barrel) average crude price, converted into local currency. Details are provided in Table 2. Current literature takes advantage of the ability of the competitive storage model to reproduce (at least partially) some of the main characteristics of commodity prices: autocorrelation, skewness, and kurtosis. However, while these characteristics are typical of international prices, they are not shared by producer prices expressed in local currencies, which are not stationary as shown in Table 3. Therefore, our estimation results presented in Table 3 are based on the difference in the log of producer prices as presented in equations ( 13) and ( 14). First, was found to be 6.176, which is close to the estimate by [START_REF] Hummels | Toward a Geography of Trade Costs[END_REF], who obtained an elasticity of approximately 5 with a standard deviation of about 2 for cereal products using a multi-sector model of trade. This estimate is also consistent with Armington elasticities proposed by McDaniel and Balistreri (2003), and classifies rice as a highly substitutable commodity. In this case, production shocks (and income changes) had rather small effects on price because of a high elasticity of substitution. Indeed, when is high, a weather shock in country has only a limited effect on producer prices in this country because consumers (in country and elsewhere) will substitute the variety of rice produced in country for varieties produced by other countries and other goods. Second, the estimate of supports nonincreasing (but close to constant) returns to scale in agriculture [START_REF] Bardhan | Size, Productivity, and Returns to Scale: An Analysis of Farm-Level Data in Indian Agriculture[END_REF]; Townsend, estimate of , the difference in column 1 was not significantly different from 0. In summary, the omission of the wage rate did not significantly affect the estimates of the structural parameters of the theoretical model. The omission of other inputs did not significantly affect the market access or exchange rate coefficients. Given that proximity to the world market is not known to be correlated with any other input prices except wages, this specification should be robust to omitted input prices.

Another problem could arise if estimation results are sensitive to the exclusion of one specific country. Table 5 contains new estimates of equation ( 14) after dropping each country in turn. No country significantly affected the estimation results. In particular, the elasticity of the market price to changes in price forecast remained high, so further inquiries must be carried out to better understand how important price expectations are in price dynamics. 7a and7b present the distribution of changes in logged prices of simulated data based on rational price expectations using actual data. The characteristics of the simulated prices were not sensitive to the choice of . However, seemed to strongly determine the shape of the distribution of simulated prices. First, all simulated series underestimated the mean price change. This clearly highlights the need for the introduction of additional mechanisms (or at least input prices) in the model, such as inventories [START_REF] Deaton | On the Behaviour of Commodity Prices[END_REF] or monetary shocks [START_REF] Frankel | Expectations and Commodity Price Dynamics: The Overshooting Model[END_REF]. Second, of the proposed values for , and generated the best characteristics in terms of volatility (represented by the coefficient of variation), skewness, and kurtosis. In order to determine which of these values for best fit the data, I generated QQ-plots (Graph 1) of simulated series against actual data for different values of , when .

While simulated series with systematically underestimated positive and negative price changes that were not close to zero, failed to predict price decreases in a satisfactory way. Even though it fails to predict the most severe price shocks, appears to be a good trade-off. In this case, producers have approximately a 50% chance of making prediction errors larger than 10% of the final price (Graph 2). This result clearly highlights the importance of price information in policies aimed at stabilizing prices.

Graph 1 -QQ-plot of simulated changes in log of prices -Rational expectations

Estimation and simulation results based on rational expectations reveal that price expectation errors are large, which suggests that policies aimed at improving the quality of price this study reinforces the current knowledge on the usefulness of the improvement of information systems through, for example, radio programs dedicated to farmers and detailing the situation of local, national and international markets [START_REF] Svensson | Getting prices right: the impact of the market information service in Uganda[END_REF]. Furthermore, since futures markets synthesize available information on future market conditions into a unique price, futures hedging could not only help peasants protect themselves against price risk, but also provide them with a new source of information, the futures price, which might help improve the quality of price expectations.

Graph 2 -Histogram of simulated expectation errors with (in proportion of the final price)

Conclusion

This article focuses on the role of price forecasting errors in a domestic price dynamics model based on international trade in varieties, rational expectations, and yield uncertainty. The theoretical model reveals that exchange rates, input prices, and price forecasts are the most influential determinants. Demand shocks and natural disasters have only a limited impact on 0 1 2 3

Density

.5 1 1.5 2 random prices because of the high elasticity of substitution between rice varieties. Empirical results support the model, confirming a high degree of substitution between rice varieties with an Armington elasticity of approximately 6 and non-increasing but close to constant returns to scale in the rice sector. Given these parameters, the model predicts that a 1% increase in price expectations results in a decrease in the market price by 1.34% and confirms the role of forecasts in price dynamics. The simulation exercise in this study was designed to determine the size of forecasting errors, and showed, assuming that forecasting errors are distributed as a log normal distribution with and , that rice producers have a 50% chance of making prediction errors larger than 10% of the final market price.

The policy implications of these findings are straightforward and twofold. First, increasing farmers access to information, such as through radio programs dedicated to providing farmers with specific information on market situations could benefit farmers by improving their access to information. Second, futures markets, which allow farmers to protect themselves against price fluctuations and provide predictions of future spot prices, may help stabilize price expectations, a necessary condition for sustainable growth.

However, the simulation exercise proposed in this paper reveals that this framework cannot fully explain the positive and negative price shocks that characterize agricultural price series. This result reveals the need for the introduction of other mechanisms in the model, such as competitive storage [START_REF] Deaton | On the Behaviour of Commodity Prices[END_REF], monetary shocks [START_REF] Frankel | Expectations and Commodity Price Dynamics: The Overshooting Model[END_REF], or trade restrictions [START_REF] Timmer | Did Speculation Affect World Rice Prices? Agricultural and Development Economics Division of the Food and Agriculture Organization of the United Nations (FAO -ESA)[END_REF], which have played an important role in recent price surges.

  log of prices Simulated change in log of price information would benefit rice producers in Asia. The problem of anticipation errors pointed in

Table 1 -

 1 Sample description

	Country	Available years Country	Available years
	Bangladesh

Table 2 -

 2 Descriptive statistics

	Variable	Mean	Standard Deviation	Minimum	Maximum
	Producer price, USD	147	51.24	30.19	384.88
	Real market potential, thousand USD	3095	5578.83	28.18	46498.73
	Fertilizer (urea) price, USD	360.35	134.37	85.86	849.72
	Agricultural GDP per capita, USD	.64	.97	.08	6.34
	Oil price, USD	17.11	9.78	1.79	35.71

Table 3 -

 3 Fisher's stationarity test for panel data (H0: unit root)

		Variable		-Statistic		P-Value
	Log of producer price, LCU	24.644		0.539
	Log of producer price, USD	54.522		0.001
	Difference in log of producer price, LCU	375.882		0.000
	Difference in log of producer price, USD	378.244		0.000
	Estimation results are presented in the first column of Table 4. As expected, all estimated
	coefficients were positive. Only exchange rates and market access coefficients were significantly
	different from zero. However, these coefficients have no economic interpretation in this
	framework due to the above-mentioned correlation with	which is part of the error term.
	Some additional steps are necessary to obtain structural parameters, and . is obtained by
	dividing	the	coefficient	of	by	the	coefficient	of
		, while requires additional calculations (with	as the coefficient of variable in
	the table of estimation results):

Table 4 -

 4 Estimation results

	Dependent variable:	(1)	(2)	(3)	(4)
		0.069**	0.071**	0.069** 0.083***
		(0.024)	(0.024)	(0 .024)	(0 .020)
		0.427*** 0.433*** 0.429*** 0.532***
		(0.044)	(0.044)	(0.045)	(0 .066)
		0.054		0.064	0.077
		(0.058)		(0 .053)	(0 .056)
		0.028	0.042		0.056
		(0.036)	(0.032)		(0 .037)
		0.146	0.160*	0.169**
		(0.085)	(0.088)	(0 .075)
	Observations	379	391	379	388
	R-squared	0.236	0.232	0.234	0.287
	(elasticity of substitution)	6.176** 6.080*** 6.217**	6.429**
		(2.039)	(1.851)	(2.058)	(2.111)
	(extent of returns to scale)	.892***	.888*** 0.892*** 0.850***
		(0.036)	(0.036)	(0.037)	(0 .028)
		-1.340*** -1.311*** -1.333*** -0.880***
		(0.240)	(0.236)	(0 .244)	(0 .232)
	Robust standard errors clustered by country in parentheses. *** p<0.01, ** p<0.05, * p<0.1

Table 5 -

 5 Robustness tests

	Country dropped			Country dropped
	Bangladesh	6.876** 0.907*** -1.426*** Nepal	6.305**	0.896*** -1.363***
		(2.749)	(0.039)	(0.253)	(2.312)	(0.039)	(0.248)
	Cambodia	6.379** 0.896*** -1.352*** Pakistan	6.009** 0.894*** -1.398***
		(2.197)	(0.036)	(0.242)	(2.141)	(0.040)	(0.265)
	China	5.553*** 0.874*** -1.260*** Philippines	5.884** 0.876*** -1.198***
		(1.724)	(0.037)	(0.238)	(1.978)	(0.038)	(0.214)
	India	5.671** 0.889*** -1.413*** Sri Lanka	7.354** 0.909*** -1.362***
		(1.867)	(0.038)	(0.237)	(3.034)	(0.038)	(0.255)
	Indonesia	7.961** 0.904*** -1.192*** Thailand	5.613*** 0.879*** -1.297***
		(2.628)	(0.039)	(0.301)	(1.742)	(0.036)	(0.243)
	Laos	6.048** 0.897*** -1.440*** Vietnam	6.180** 0.894*** -1.368***
		(2.073)	(0.038)	(0.258)	(2.107)	(0.037)	(0.243)
	Malaysia	5.536*** 0.880*** -1.328***
		(1.604)	(0.036)	(0.251)
	Robust standard errors clustered by country in parentheses. *** p<0.01, ** p<0.05, * p<0.1

Table 6 -

 6 Annual rice production and losses arising as a consequence of natural disasters in thePhilippines between 1991Philippines between -2000

			Typhoons and Floods	Drought	Total losses
	Year	Total Palay production (thd. tons)	Loss in production (thd. tons)	% Loss of production	Loss in production (thd. tons)	% Loss of production	Loss in production (thd. tons)	% Loss of production
	1991	9673	152	1.6	16	0.2	168	1.7
	1992	9129	11	0.1	79	0.9	90	1
	1993	9434	246	2.6	-	0	246	2.6
	1994	10538	112	1.1	48	0.5	160	1.5
	1995	10541	328	3.1	45	0.4	373	3.5
	1996	11284	73	0.6	1	0	74	0.7
	1997	11269	75	0.7	15	0.1	90	0.8
	1998	8555	1048	12.3	462	5.4	1511	17.7
	1999	11787	321	2.7	1	0	322	2.7
	2000	12389	390	3.1	-	0	390	3.1
	Total	104599	2757	2.6	667	1.5	3424	3.3
	Source: (Garcia-Rincon and Virtucio 2008)
		Tables

Table 7a -

 7a Summary statistics of observed and simulated price changes

		Observed data	Simulated data
	Mean	0.070	0.055
	Standard deviation	0.161	0.065
	Coefficient of variation	2.285	1.162
	Skewness	0.757	2.907
	Kurtosis	6.791	19.070

Table 7b -

 7b Summary statistics of simulated price changes

	Mean	0. 055	0.055	0.056	0.056	0.058
	Std. dev.	0. 065	0.077	0.103	0.142	0.176
	Coef. of var.	1.174	1.388	1.830	2.527	3.042
	Skewness	2.857	2.178	0.970	-0.002	0.220
	Kurtosis	18.822	14.596	6.386	3.385	4.032
	Mean	0.055	0.055	0.054	0.054	0.057
	Std. dev.	0.065	0. 076	0.099	0.121	0.169
	Coef. of var.	1.172	1.370	1.820	2.230	2.972
	Skewness	2.807	1.544	0.746	0.497	0.335
	Kurtosis	18.100	8.989	6.031	4.573	3.563
	Mean	0.056	0.056	0.056	0.056	0.055
	Std. dev.	0.065	0.077	0.103	0.135	0.180
	Coef. of var.	1.177	1.382	1.856	2.397	3.263
	Skewness	2.736	1.811	1.087	-0.008	-0.195
	Kurtosis	17.960	10.961	8.108	3.953	3.166

The study is narrowed to the 1965-2003 period due to the lack of availability of market access data prior to and after this interval(Head and Mayer

2011).

 Kirsten, and Vink 1998). Given these parameters, I estimated the elasticity of the final price with respect to price expectations as . I obtained a surprisingly large effect of price expectations on market price of .

Hence, a 1% increase in the forecasted price resulted in a fall of the market price by 1.34%, all things being equal. This result clearly highlights the importance of policies aimed at improving the quality of price information and stabilizing price expectations. Apart from these results, none of the included input prices appeared to have a significant effect, but they all had positive coefficients, in line with the predictions of the theoretical model.

Before further analyzing the various implications of these results, I used several methods to test their robustness. First, the scarcity of input price data was a major concern for the proposed estimation strategy. The results in columns 2-4 in Table 4 were obtained by dropping input prices in order to verify the stability of the coefficients of exchange rates and market potential variables, which are of crucial importance in the calculations of the structural parameters and and the elasticity of market price to price forecasts. If the estimated coefficients are sensitive to the inclusion or exclusion of one input price, then the omission of other input prices should bias the results. Only the exclusion of wage rate data significantly affected the coefficients of interest. This is by no means a surprise given the growing literature on the link between market access and wages. Many theoretical and empirical analyses have found a positive link between proximity to world markets and wages, which was confirmed by comparing the differences in coefficients between column 4 and all the other columns in Table 4 since the effect of increases when the wage rate drops. This change in the coefficient of is counterbalanced by an increase in the coefficient of , so that neither nor were significantly affected by the omission of the wage rate. Despite a decrease in the point

The present model does not allow us to decompose the error term between weather shocks and forecasting errors. However, equations ( 13) and ( 14) can give some information on the correlation between and .

First, according to equation ( 13), . Second, according to equation ( 14), . We obtain:

(15)

Reformulating this last equation gives:

(16)

Given this correlation, the objective was to simulate forecasting errors correlated with production shocks assuming different standard deviations for forecasting errors and find which provided the best fit to the data. To do so, I simulated both weather shocks and forecasting errors, assuming that and and tested for different values for and . To do so, I chose and in accordance with the data presented in


Mridul Aanjaneya
email: aanjneya@stanford.edu

Frederic Chazal
email: frederic.chazal@inria.fr

Daniel Chen
email: danielc@cs.stanford.edu

Marc Glisse
email: marc.glisse@inria.fr

Leonidas Guibas
email: guibas@cs.stanford.edu

Dmitriy Morozov
email: dmitriy@mrzv.org

Metric Graph Reconstruction from Noisy Data

. Building on the notions of correspondence and Gromov-Hausdorff distance in metric geometry, we describe a model for such data sets as an approximation of an underlying metric graph. We present a novel algorithm that takes as an input such a data set, and outputs the underlying metric graph with guarantees. We also implement the algorithm, and evaluate its performance on a variety of real world data sets.

Introduction

Motivation. Large-scale geometric data sets are becoming widely available, whether from high-bandwidth sensors or from massive simulations of physical processes. All across science, engineering, medicine, and defense, there is a real need to analyze, understand, and extract useful information out of such massive geometric data. Much of this data is noisy, contains outliers, has missing parts, and does not have a manifold structure or even a consistent dimension -raising many difficult statistical, geometric, and algorithmic problems in its analysis. In this paper, we focus on a simple, but important setting of mixeddimension geometric data, namely a setting where the underlying space of the data can be viewed of as a metric graph [START_REF] Kuchment | Quantum graphs I. Some basic structures[END_REF], which is a 1-D stratified space consisting of just 0-D strata (vertices) and 1-D linear strata (edges or loops), glued together in some fashion, see Figure 1(a).

Branching filamentary structures, which can be naturally viewed of as metric graphs, appear in a wide variety of real-world data sets, both in settings where the data arises embedded in Euclidean space, as well as in situations where the host space is less intuitive and only local metric information may be available. For example large-scale collections of GPS traces for vehicles or pedestrians are becoming widely available (see e.g., [2]) and can be used to provide a variety of location-aware services. Their movement patterns tend to follow a branching structure which can be modeled as a metric graph. Earthquake faults are intimately connected with plate tectonics and tend to follow filamentary structures as they arise along the boundaries of such plates (see e.g., [START_REF]Earthquake search[END_REF]). In nuclear physics, high-energy particles move along filamentary trajectories and there is often the need to track their motion [START_REF] Abramowicz | Orientation selective neural network for cosmic muon identification[END_REF]. In materials science, stresses can cause material cracks that propagate along branching structures formed by linear paths; their detection is an important research problem [START_REF] Qaddoumi | Microwave detection of stress-induced fatigue cracks in steel and potential for crack opening determination[END_REF]. Many defense applications require the extraction of road networks from synthetic aperture radar (SAR) images [START_REF] Tupin | Detection of linear features in SAR images: Application to road network extraction[END_REF]. In astronomy, filamentary structures in galaxies are of great interest (e.g., [START_REF] Choi | Tracing the filamentary structure of the galaxy distribution at z ∼ 0.8[END_REF]) for cosmological studies. This is not to mention networks formed by blood vessels in the body for anatomy, river systems in geography, and many other examples.

Branching structures are also quite common in more abstract settings, though sometimes one has to look at such data with a coarser lens before it becomes apparent. For instance, communication networks can be regarded as large graphs in which certain dominant pathways define the major arteries connecting network hubs. Recently, Heath et. al. [START_REF] Heath | Image webs: Computing and exploiting connectivity in image collections[END_REF] built large graphs from image collections, by linking together images with partial shared content. Data sets of interest here include collections of images acquired by a mobile agent along its path, as in Google Streetview. In such cases, at a coarse scale, the connectivity among the images reflects the mobility of the capturing agent, naturally giving rise to branching filamentary structures. Extraction of this underlying structure can provide a useful map for understanding the image data, navigating through it, or for answering certain queries.

Reconstruction Problem. While there has been a great deal of prior work on both topological and geometric reconstruction of geometric data sets under varying sampling conditions, our emphasis is on an intermediate level of reconstruction, what we term metric reconstruction -a largely unexplored domain. The input to our algorithm is a metric space (Y, d Y ) that is close to a much simpler metric graph (X, d X ) in a sense that we make precise in the Section 3. (Y, d Y ) can be constructed from raw data in various ways: in some cases, we construct a neighborhood graph on the raw data, and use the shortest path as the distance; in other cases, the metric is given to us directly. Note that this implies that our reconstruction is aimed at capturing the intrinsic structure of the data and is somewhat oblivious to its extrinsic embedding, wherever that is available. Our goal is then to extract a metric graph ( X, d X) that has the same topology as (X, d X ), and a map φ : Y → X that approximately preserves distances.

Experiments. In addition to theoretical reconstruction results with performance guarantees, we study experimentally the performance of our algorithm on a variety of data sets from different applications, including data in which an embedding is given (GPS traces, earthquake data, astronomical data), as well as data in which only metric information is available (Image Webs). In all these cases our compact metric approximation provides a much more manageable representation of the structure of the original data -far easier to visualize, navigate, and manipulate than the original. Our metric guarantees allow us to further exploit this representation by running graph algorithms in this compact representation in lieu of the original graph. As an example, we used the compressed graph to perform shortest path queries, resulting in significant speedups on some data sets. Related Work. Our work is related to contributions by several different communities. On the one side, the statistics community has investigated the problem of extracting filamentary structures from point cloud data, starting with the seminal work of Arias-Castro et al. [START_REF] Arias-Castro | Adaptive multiscale detection of filamentary structures in a background of uniform random points[END_REF] based on counting membership in multiscale anisotropic strips. Subsequent approaches exploit gradient descent or medial axis ideas [START_REF] Genovese | On the path density of a gradient field[END_REF][START_REF] Genovese | Nonparametric Filament Estimation[END_REF]. All these, however, aim mostly at the extraction of isolated filaments, focus on how to deal with outlier data, and do not pay serious attention to the global branching structures the filaments form. Also, they all assume an extrinsic embedding of the data. On the other side, there has been extensive work in the computational geometry community on curve reconstruction, which is the problem of computing a polygonal curve that approximates well a curve sampled by a given point set -several algorithms have been proposed for this problem [START_REF] Amenta | The crust and the β-skeleton: Combinatorial curve reconstruction[END_REF][START_REF] Dey | Curve reconstruction: Connecting dots with good reason[END_REF][START_REF] Dey | Reconstructing curves with sharp corners[END_REF]. Unfortunately, it is hard to extend these methods to our setting, since they also view 0-dimensional strata, which exhibit non-manifold behavior, as singularities and try to avoid them as much as possible. While geometric reconstruction is not our goal, as in that work, we aim to be able to prove certain quality guarantees on the metric reconstruction we attain, under appropriate sampling conditions. Finally, Chen et al. [START_REF] Chen | Road network reconstruction for organizing paths[END_REF] recently considered a related problem of reconstructing a road network from a given collection of path traces. They designed an algorithm with guarantees without making heavy assumptions on the distribution of input paths. However, the assumptions they use are stronger than desired in many practical applications. In particular, their method depends on an embedding of the data and sequential path information.

We end by remarking that dimension reduction has been a topic of much study in the machine learning and data analysis communities. When data is given in parametric form, i.e., as points in a (possibly high dimensional) Euclidean space and the goal of dimension reduction is distance preservation, many wellknown methods exist based on random projections as suggested by the Johnson-Lindenstrauss lemma, or by locality sensitive hashing (LSH). This paper addresses "dimension reduction" for distance preservation in the case where the metric is given by the shortest path distance on a large but special type of graphone that contains few but large linear structures. As we show, this type of metric reconstruction raises interesting new mathematical problems and is applicable to many types of geometric data.

Preliminaries

Recall that a metric space is a pair (X, d) where X is a set and d : X × X → R + is a symmetric function satisfying (1): d(x, x ′ ) = 0 if and only if x = x ′ and (2): d(x, x ′′ ) ≤ d(x, x ′ ) + d(x ′ , x ′′ ). Two spaces (X, d X ) and (Y, d Y ) are isometric if there exists a bijection φ : X → Y that preserves the distances, namely: d Y (φ(x), φ(x ′ )) = d X (x, x ′ ) for all x, x ′ ∈ X. The space of isometry classes of metric spaces is endowed with the Gromov-Hausdorff distance [START_REF] Gromov | Metric Structures for Riemannian and Non-Riemannian Spaces[END_REF] whose definition can be given using the notion of ε-correspondences ([6] Thms 7.3.25 and 7.3.30).

Definition 1. A correspondence between (X, d X ) and (Y, d Y ) is a set C ⊂ X × Y
such that for any x ∈ X (resp. y ∈ Y), there exists y ∈ Y (resp. x ∈ X) such that (x, y) ∈ C. When x, y are such that (x, y) ∈ C, we say that x and y are paired in C. Given ε > 0, C is an ε-correspondence if for any (x, y),

(x ′ , y ′ ) ∈ C, |d X (x, x ′ ) -d Y (y, y ′ )| ≤ ε. The Gromov-Hausdorff distance d GH (X, Y)
is the infimum of the ε ≥ 0 such that there exists an ε-correspondence between (X, d X ) and (Y, d Y ).

An ε-correspondence between X and Y can be seen as an ε-approximation of X by Y (and reciprocally). However, in many applications, data only comes with locally correct approximate metric information. For example, for a data set sampling a road network the Euclidean distance between data points provides a suitable approximation of the metric of the underlying network only locally. So in this paper, we use a more local and weaker notion of correspondence: given positive numbers ε, R, we say that (Y,

d Y ) is an (ε, R)-approximation of a metric space (X, d X ) if there exists a correspondence C ⊂ X × Y such that (x, y), (x ′ , y ′ ) ∈ C, min(d X (x, x ′ ), d Y (y, y ′ )) ≤ R =⇒ |d X (x, x ′ ) -d Y (y, y ′ )| ≤ ε
Notice that this latter notion is strictly weaker than the notion of ε-correspondence. In particular, the existence of an (ε, R)-correspondence between X and Y does not bound d GH (X, Y), as shown in the following example: let X ⊂ R 2 be the half-circle {x 2 + y 2 = 1, y ≤ 0} endowed with the geodesic distance and let Y = X be endowed with the restriction of the Euclidean distance. For any ε > 0, the diagonal C = {(x, x) : x ∈ X} ⊂ X × Y is an (ε, O(ε1/3 ))-correspondence, but the diameters of X and Y are respectively equal to π and 2, showing that d GH (X, Y) ≥ π -2 > 0. Nevertheless, (ε, R)-approximations give rise to global approximations with respect to d GH when the approximated space is a path metric space, defined as follows: Definition 2. A metric space (X, d X ) is a path metric space if the distance between any pair of points is equal to the infimum of the lengths of the continuous curves joining them 1 . Equivalently (X, d X ) is a path metric space if and only if for any x, y ∈ X and any ε > 0 there exists z ∈ X such that max(d X (x, z), d X (y, z)) ≤ 1 2 d X (x, y) + ε [START_REF] Gromov | Metric Structures for Riemannian and Non-Riemannian Spaces[END_REF]. Then, we can obtain the following bound on the Gromov-Hausdorff distance: Lemma 3. Let (X, d X ) be a path metric space, (Y, d Y ) an (ε, R)-approximation of X and assume that Y has the following property:

(⋆) for any y, y ′ ∈ Y there exists a sequence y 0 = y, y 1 , . . . y n-1 , y n = y ′ such that for all i = 0, . .

. , n -1, d Y (y i , y i+1 ) ≤ R and d Y (y, y ′ ) = n-1 i=0 d Y (y i , y i+1 ). If C ⊂ X × Y is an (ε, R)-correspondence, then for any (x, y), (x ′ , y ′ ) ∈ C we have |d Y (y, y ′ ) -d X (x, x ′ )| ≤ min(d X (x, x ′ ), d Y (y, y ′ )) R/2 + 1 ε. In particular, d GH ((X, d X ), (Y, d Y )) ≤ diam(X) R/2 + 1 ε where diam(X) is the diameter of X.
Proof. (x, y) and (x ′ , y ′ ) ∈ C are given. By hypothesis, there exists a sequence y 0 = y, y 1 , . . . y n-1 , y n = y ′ such that for all i = 0, . .

. , n -1, d Y (y i , y i+1 ) ≤ R and d Y (y, y ′ ) = n-1 i=0 d Y (y i , y i+1 ). As a first remark, notice that in this sequence, if i < j, d Y (y i , y j ) = j-1 k=i d Y (y k , y k+1 )
Indeed, using the triangle inequality:

d Y (y, y ′ ) ≤ i-1 k=0 d Y (y k , y k+1 ) + d Y (y i , y j ) + n-1 k=j d Y (y k , y k+1 ) ≤ n-1 k=0 d Y (y k , y k+1 ) = d Y (y, y ′ ).
In property (⋆), we can further assume that d Y (y i , y i+2 ) > R. If d Y (y i , y i+2 ) ≤ R, we can remove y i+1 from the sequence, and the previous remark shows that the properties are still satisfied. In particular, this implies that d Y (y, y ′ ) > n-1 2 R. Because each y i has a corresponding (non-unique) x i ∈ X in C, we have:

d X (x, x ′ ) ≤ n-1 i=0 d X (x k , x k+1 ) ≤ n-1 i=0 d Y (y k , y k+1 ) + nε ≤ d Y (y, y ′ ) + d Y (y, y ′ ) R/2 + 1 ε.
A simple computation shows that this implies:

d X (x, x ′ ) < d Y (y, y ′ ) + d X (x, x ′ ) R/2 + 1 ε.
Now X almost satisfies (⋆). Indeed, by recursively splitting the intervals of length more that R, for any ε ′ > 0, we construct a sequence

x 0 = x, x 1 , . . . , x n = x ′ such that d X (x i , x i+1 ) ≤ R, d X (x i , x i+2 ) ≥ R -ε ′ and n-1 i=0 d X (x i , x i+1 ) ≤ d Y (y, y ′ ) + ε ′ .
We derive as before:

d Y (y, y ′ ) ≤ d X (x, x ′ ) + d X (x, x ′ ) + ε ′ R/2 -ε ′ + 1 ε
and since it is true for all ε ′ > 0:

d Y (y, y ′ ) ≤ d X (x, x ′ ) + d X (x, x ′ ) R/2 + 1 ε
which again implies:

d Y (y, y ′ ) < d X (x, x ′ ) + d Y (y, y ′ ) R/2 + 1 ε
and finally:

|d Y (y, y ′ ) -d X (x, x ′ )| ≤ min(d X (x, x ′ ), d Y (y, y ′ )) R/2 + 1 ε
In this paper, we assume that our input is an (ε, R)-approximation of a specific type of path metric space, known as a metric graph [START_REF] Kuchment | Quantum graphs I. Some basic structures[END_REF]:

Definition 4. A metric graph is a path metric space (X, d X ) that is homeomorphic to a 1-dimensional stratified space (see Figure 1(a)). A vertex of X is a 0-dimensional stratum of X and an edge of X is a 1-dimensional stratum of X 2 .
It is useful to note that edges are isometric to finite length intervals in the real line.

Problem Definition

Let (Y, d Y ) be an (ε, R)-approximation of a metric graph (X, d X ) that has a shortest edge length of b. Without loss of generality, we will assume that X is connected. Note that our definition of (ε, R)approximation is essentially a worst-case noise model for the data that does not rely on further distributional assumptions.

(ε, R)-approximation of (X, d X ). Our goal is to design an algorithm to reconstruct from (Y, d Y ) a space ( X, d X) that is homeomorphic to (X, d X )
. Furthermore, we define distances on ( X, d X) that approximate those of (X, d X ) and return a map φ : Y → X that approximately preserves distances. Although we frame this objective as a reconstruction problem, in practice, our algorithm can be used to find a much simpler metric graph ( X, d X) approximating the input space (Y, d Y ), and achieving guarantees when (Y, d Y ) is an approximation of a suitable metric graph (X, d X ).

Overview of Algorithm and Guarantees

In addition to the input metric (Y, d Y ), which, as mentioned previously, is an (ε, R)-approximation of an underlying metric graph (X, d X ), our algorithm also takes a parameter r that roughly corresponds to the scale at which we look at the data. For noisier data, we would generally use a larger r, while to capture smaller features, we would choose a smaller r. Our analysis will exhibit a range of values for r that result in a correct reconstruction depending on both ε and R, as well as on b, the length of the shortest edge in X. In practice, we do not know these values, but our implementation always outputs a suitable metric graph ( X, d X) for which we can check the distortion of the metric from (Y, d Y ). Hence, we are able to try values of r until we obtain a suitable and simple approximation of (Y, d Y ).

Recall that there is an (ε, R)-correspondence between our input metric (Y, d Y ) and its underlying metric graph (X, d X ). The algorithm proceeds in two steps. First it begins by labeling as "branch points" the points of Y paired under this correspondence to a point in X that is close to a vertex and labeling the rest of the points of Y as "edge points". Then, the algorithm uses these labels to reconstruct a new metric graph X that is homeomorphic to X and estimates distance preserving maps from Y to X. For ease of reference, the pseudocode of our algorithm is given in Algorithm 1.

The following results show that if (Y, d Y ) is a sufficiently good approximation of (X, d X ) then the reconstructed graph ( X, d X) is homeomorphic and almost isometric to (X, d X ).

Theorem 5 (Topological Reconstruction).

If the length b of the shortest edge of X is larger than 16r and 15ε/2 < r < min(R/4, 3(b -2ε)/5) then the reconstructed graph X is homeomorphic to X.

Theorem 6 (Metric Reconstruction).

Under the assumptions of Theorem 5 there exists a homeomorphism φ : X → X such that for any

Algorithm 1 Metric Graph Reconstruction

Require: Metric space (Y, d Y ) approximating metric graph (X, d X ) and parameter r > 0. Ensure: Metric graph ( X, d X)

1: Labeling points as edge or branch 2: for all y ∈ Y do if degr(y) = 2 then 6:

Label y as a edge point. Label y as a preliminary branch point.

9:

end if 10: end for 11: Label all points within distance 2r from a preliminary branch point as branch points.

12: Let E be the points of Y labeled as edge points. 13: Let V be the points of Y labeled as branch points. 14: Reconstructing the Graph Structure 15: Compute the connected components of the Rips-Vietoris graphs R2r(E) and R2r(V). 16: Let the connected components of R2r(V) be the vertices of the reconstructed graph X. 17: Let there be an edge between vertices of X if their corresponding connected components in R2r(V) contain points at distance less than 2r from the same connected component of R2r(E). 18: Reconstructing the Metric 19: To each edge ê of X assign a length equal to the diameter of the corresponding connected component of R2r(E) plus 4r.

x,

x ′ ∈ X, (1-κ)d X (x, x ′ ) ≤ d X(φ(x), φ(x ′ )) ≤ (1+κ ′ )d X (x, x ′ ) with κ = 10r 3b +( 5 b + 2 R )ε and κ ′ = ( 3 b + 2 R )ε.
The proofs of these results, as well as a more detailed discussion of the algorithm, are given in the next section where an easy to compute map with low-metric distortion between (Y, d Y ) and ( X, d X) is also provided.

Analysis and Proofs

In this section we assume that the assumptions of Theorems 5 and 6 are satisfied.

Labeling points as edge or branch

First notice that the classification of a point x ∈ X as a vertex or a point on an edge is determined by the number of connected components of a small intrinsic sphere centered at x (see Figure 1(b)). To label a point y ∈ Y as either a branch point or an edge point, our algorithm considers the intrinsic spherical shells B Y (y, 5r/3) \ B Y (y, r) around y and constructs a Rips-Vietoris graph with parameter 4r/3 on the points of Y inside the spherical shell. Then, it records the number of connected components of this graph as the r-degree deg r (y): r) and edges connecting all the pairs of vertices at distance less than 4r/3 from each other.

Definition 7. Let (Y, d Y ) be an (ε, R)-approximation of X. Given 0 < r < R/2, the r-degree deg r (y) of a point y ∈ Y is the number of connected components of the Rips-Vietoris graph R 4r/3 (B Y (y, 5r/3) \ B Y (y, r)) with vertex set B Y (y, 5r/3) \ B Y (y,
Intuitively, it is easy to imagine that if deg r (y) = 2, then y corresponds to a point on X close to a vertex, whereas if deg r (y) = 2, y corresponds to a point on X far from a vertex.

Theorem 8 (Degree Inference Theorem). Let (Y, d Y ) be an (ε, R)-approximation of X. Let C ⊂ X × Y be an (ε, R)-correspondence between X and Y, let (x, y) ∈ C. i) If the distance d 0 from x to any vertex of X is larger than 17 2 ε, then for 9 2 ε < r < min( R 2 , 3(d0-ε) 5
), deg r (y) is equal to the degree of x in X (i.e. 2). Moreover the pairwise distances between the connected components of the Rips-Vietoris graph are lower bounded by 2r -3ε. ii) If x is at distance less than ε from a vertex x 0 of X and if the length l 0 of the shortest edge adjacent to x 0 is larger than 27 2 ε then for 15 2 ε < r < min( R 2 , 3(l0-2ε)

5

), deg r (y) is equal to the degree of x 0 in X. Moreover the pairwise distances between the connected components of the Rips-Vietoris graph are lower bounded by 2r -5ε.

Proof. This theorem is a consequence of Theorem 10 with α set to 2/3.

Our proof of the theorem follows from a more general theorem that considers a variation of the rdegree dependent on an extra parameter 0 ≤ α ≤ 1. This variation allows us to adjust the radius of the ball B Y (y, 5r/3) to any (1 + α)r for different guarantees:

Let (Y, d Y ) be an (ε, R)-approximation of X. Let C ⊂ X × Y be an (ε, R)-correspondence between X and Y , let (x, y) ∈ C and let 0 < α < 1. i) If the distance d 0 from x to any vertex of X is larger than (3 max( 1+α α , 1+α 2(1-α) ) + 1)ε then for 3 max( 1 α , 1 2(1-α) )ε < r < min( R 2 , d0-ε 1+α ), deg r,α (
y) is equal to the degree of x in X (i.e. 2). Moreover the pairwise distances between the connected components of the Rips-Vietoris graph are lower bounded by 2r -3ε. ii) If x is at distance less than ε from a vertex x 0 of X and if the length l 0 of the shortest edge adjacent to x 0 is larger than

(max( 3(1+α) α , 5(1+α) 2(1-α) ) + 1)ε then for max( 3 α , 5 2(1-α) )ε < r < min( R 2 , l0-2ε 1+α ), d r,α (y
) is equal to the degree of x 0 in X. Moreover the pairwise distances between the connected components of the Rips-Vietoris graph are lower bounded by 2r -5ε.

Proof. First remark that if (x ′ , y ′ ) ∈ C is such that y ′ ∈ B Y (y, (1 + α)r) \ B Y (y, r) then x ′ ∈ B X (x, (1 + α)r + ε) \ B X (x, r -ε). i) Since r > ε and (1 + α)r + ε < d 0 , B X (x, (1 + α)r + ε)\ B X (x, r -ε) is included in the edge containing
x and has exactly 2 connected components. Moreover, these two connected components are at distance 2(r -ε). Now, if (x ′ , y ′ ), (x ′′ , y ′′ ) ∈ C are such that y ′ , y ′′ ∈ B Y (y, (1 + α)r) \ B Y (y, r) and d Y (y ′ , y ′′ ) < 2αr then d X (x ′ , x ′′ ) < 2αr + ε and, since r > 3ε 2(1-α) , it follows that x ′ and x ′′ are in the same connected component of B X (x, (1 + α)r + ε) \ B X (x, r -ε).

Reciprocally, if (x ′ , y ′ ), (x ′′ , y ′′ ) ∈ C are such that x ′ , x ′′ are in the same connected component of B X (x, (1 + α)r + ε) \ B X (x, r -ε), then d X (x ′ , x ′′ ) ≤ αr + 2ε and d Y (y ′ , y ′′ ) ≤ αr + 3ε < 2αr since αr > 3ε.

As a consequence, the Rips-Vietoris graph with parameter 2αr and vertex set B Y (y, (1+α)r)\B Y (y, r) has at most two connected components. To prove that it has exactly two connected components one just needs to check that each connected component

K of B X (x, (1 + α)r + ε) \ B X (x, r -ε) contains a point x ′ such that there exists y ′ ∈ B Y (y, (1 + α)r) \ B Y (y, r) satisfying (x ′ , y ′ ) ∈ C: let x ′ be the point of K such that d X (x, x ′ ) = (1+α/2)r and let (x ′ , y ′ ) ∈ C. Then, since αr > 2ε, d Y (y, y ′ ) ≤ (1+α/2)r +ε < (1+α)r and d Y (y, y ′ ) ≥ (1 + α/2)r -ε > r.
ii) This is almost the same proof as for i) except that since x is not a vertex, but at distance at most ε from a vertex we have to slightly change the constraint on r.

Let y 0 ∈ Y be such that (x 0 , y 0 ) ∈ C. From the definition of (ε, 0)-approximation we have d Y (y, y 0 ) < 2ε.

Since r > 2ε and (1 + α)r + 2ε < l 0 , B X (x, (1 + α)r + ε) \ B X (x, r -ε) has exactly d connected components, each included in different edges adjacent to x 0 , where d is the degree of x 0 . Moreover these connected components are at distance a least 2(r -2ε) from each other. Now, if (x ′ , y ′ ), (x ′′ , y ′′ ) ∈ C are such that y ′ , y ′′ ∈ B Y (y, (1 + α)r) \ B Y (y, r) and d Y (y ′ , y ′′ ) < 2αr then d X (x ′ , x ′′ ) < 2αr + ε and, since r > 5ε 2(1-α) , it follows from claim 1 that x ′ and x ′′ are in the same connected component of B X (x, (1 + α)r + ε) \ B X (x, r -ε).

Reciprocally, if (x ′ , y ′ ), (x ′′ , y ′′ ) ∈ C are such that x ′ , x ′′ are in the same connected component of B X (x, (1 + α)r + ε) \ B X (x, r -ε), then d X (x ′ , x ′′ ) ≤ αr + 2ε and d Y (y ′ , y ′′ ) ≤ αr + 3ε < 2αr since αr > 3ε.

As a consequence, the Rips-Vietoris graph with parameter 2αr and vertex set B Y (y, (1+α)r)\B Y (y, r) has at most d connected components. To prove that it has exactly d connected components one just needs to check that each connected component K of B X (x, (1 + α)r + ε) \ B X (x, r -ε) contains a point x ′ such that there exists y ′ ∈ B Y (y, (1 + α)r) \ B Y (y, r) satisfying (x ′ , y ′ ) ∈ C: let x ′ be the point of K such that d X (x, x ′ ) = (1+α/2)r and let (x ′ , y ′ ) ∈ C. Then, since αr > 2ε, d Y (y, y ′ ) ≤ (1+α/2)r +ε < (1+α)r and d Y (y, y ′ ) ≥ (1 + α/2)r -ε > r.

This result motivates the choice of the value α = 2/3 in Theorem 8: this is the value that minimizing the bound (3 max( 1+α α , 1+α 2(1-α) ) + 1)ε in i) that controls the size of the expansion procedure in Section 5.2.

Reconstructing the graph structure

We now describe the reconstruction procedure. Given 15ε/2 < r < min(R/4, 3(b -2ε)/5), recall that we first label the points y ∈ Y as branch or edge depending on deg r (y): y is labelled as an edge point if deg r (y) = 2, and labelled as a branch point otherwise. The following result is an immediate consequence of Theorem 8.

Lemma 11. If y ∈ Y is paired in C to a point x at distance at most ε from a vertex of X then y is labeled as a branch point by the procedure above. If y ∈ Y is paired in C to a point x at distance at least 5r/3 + ε from any vertex of X then y is labeled as an edge point.

The points of Y paired to points in X that are at distance between ε and 5r/3 + ε from a vertex of X can be "incorrectly" labeled as branch points. It is not possible to distinguish these fuzzy points from the data Y only, so we force them to be branch points using the following expansion procedure: all points y ∈ Y that are at distance at most 2r from a point labeled as branch are promoted to branch points.

To prove that after this expansion all the fuzzy points are labeled as branch, notice that if y ∈ Y is now labeled as an edge then it is at distance at least 2r from any point y ′ ∈ Y labeled as branch before the expansion procedure. It follows that for any pair (x, y) ∈ C, x is at distance more than 2r -ε > 5r/3 + ε (since r > 15ε/2) from a vertex of X.

Corollary 12. Let (x, y) be a pair in C. If x is at distance at least 11r/3 + 2ε from any vertex of X, then after the expansion procedure, y is labeled as an edge. Reciprocally, if y is labeled as an edge after the expansion procedure, then x is at distance at least 2r -ε from a vertex of X. Now to recover the connectivity of X, we group the branch points (resp. the edge points) in clusters, each corresponding to a vertex (resp. an edge) of X. For that, we consider the Rips-Vietoris graph R 2r (V) (resp. R 2r (E)) of parameter 2r built on top of the set V ⊂ Y of branch points (resp. the set E ⊂ Y of edge points). Lemma 13. If the length b of the shortest edge of X is larger than 16r then the connected components of R 2r (V) are in one-to-one correspondence with the vertices of X and the connected components of R 2r (E) are in one-to-one correspondence with the edges of X.

Proof. If y ∈ Y is a branch point and (x, y) ∈ C, then there exist (x ′ , y ′ ), (x 0 , y 0 ) ∈ X such that x 0 is a vertex in X, d X (x 0 , x ′ ) ≤ 5r/3 + ε and d X (y ′ , y) ≤ 2r. It follows that d Y (y 0 , y ′ ) ≤ 5r/3 + 2ε ≤ 2r. So y and y 0 are in the same connected component of R 2r (V) and d Y (y, y 0 ) ≤ 4r. As a consequence, any connected component of R 2r (V) contains at least one point paired with a vertex of X.

Now if (x 1 , y 1 ) ∈ C is such that x 1 is another vertex of X, then Lemma 3 implies that d Y (y 0 , y 1 ) ≥ d X (x 0 , x 1 ) -( 2d X (x 0 , x 1 ) R + 1)ε ≥ 4 5 b
where to get the last inequality we used that d X (x 0 , x 1 ) ≥ b, R > 15ε and b > 15ε. Assume that y 0 and y 1 are in the same connected component of R 2r (V). Then there exists a path joining y 0 to y 1 in this component and since x 0 = x 1 , there exists a branch point y ′ ∈ Y along this path such that b/2 -r ≤ d Y (y ′ , y 0 ) ≤ b/2 + r. According to Lemma 3 for any x ′ ∈ X such that (x ′ , y ′ ) ∈ C, we have

d Y (y ′ , y 0 )(1 - 2ε R ) -ε ≤ d X (x ′ , x 0 ) ≤ d Y (y ′ , y 0 )(1 + 2ε R ) + ε
Using again that ε/R < 1/15 and ε < b/15 we get 11 30 b -13 15 r ≤ d X (x ′ , x 0 ) ≤ 19 30 b + 17 15 r and since b is the length of the shortest edge of X, the distance between x ′ and any vertex of X is at least min( r. Since b > 16r, one deduces from the corollary 12 that x ′ is an edge point: a contradiction. As a consequence, the points of any connected component of R 2r (V) can be paired with at most one vertex of X. This proves that the connected components of R 2r (V) are in one-to-one correspondence with the vertices of X.

To prove the second part of the lemma, first notice that since b > 16r for any edge of X there exists a point at distance at least 8r from any vertex of X. As a consequence, any y ∈ Y such that (x, y) ∈ C is labeled as an edge point showing that E contains points from all the edges of X. Now if (x, y), (x ′ , y ′ ) ∈ C are such that y, y ′ ∈ E and x, x ′ are not in the same edge of X, then any shortest path joining x to x ′ has to meet a vertex x ′′ of X. So for any sequence (x 0 , y 0 ) = (x, y), (x 1 , y 1 ),

• • • (x n , y n ) = (x ′ , y ′ ) ∈ C such that y 0 = y, y 1 . . . , y n = y ′ is joining y to y ′ in R 2r (Y) there exists i ∈ {1, • • • , n -1} such that d X (x ′′ , x i ) ≤ 1 2 (2r + ε).
It follows that y i is a branch point and y and y ′ cannot be in the same connected component of R 2r (E). Reciprocally, if x, x ′ are in the same edge e of X, they both are at distance at least 2r -ε from the end points of e and from any point paired to a point labeled as branch before the expansion procedure. So, if (x ′′ , y ′′ ) ∈ C is such that x ′′ ∈ e is contained in the interval defined by x and x ′ and is at distance larger than 2r -ε from x and x ′ , then the distance from x ′′ to any point paired to a branch point before the expansion procedure is at least 4r -2ε > 11r/3 + 2ε (since r > 15ε/2). Therefore, y ′′ is an edge point. As a consequence, there exists a sequence y 0 = y, y 1 . . . , y n = y ′ of edge points that are all paired to points in the edge e such that d Y (y i , y i+1 ) ≤ 2r for i = 0, . . . n -1, proving that y and y ′ are in the same connected component of R 2r (E). It follows that the connected components of R 2r (E) are in one-to-one correspondence with the edges of X. Now recall that X is built as follows: we create a vertex for each connected component of R 2r (V); we create an edge between two vertices if each of the two corresponding components contains at least one point at distance less than 2r from the same connected components of R 2r (E). From Lemma 13 we then deduce the Topological Reconstruction Theorem 5.

Reconstructing the metric

We begin with the proof of Theorem 6:

Proof of Theorem 6. The proof consists of showing the existence of a (1 + κ ′ )-Lipschitz homeomorphism φ : X → X with inverse (1 -κ) -1 -Lipschitz. To this end, we proceed with each edge separately. Let ê be an edge of X, and let y 0 , y 1 be two points in the corresponding connected component in R 2r (E) such that d Y (y 0 , y 1 ) is equal to the diameter of this component. Denoting e the edge of X corresponding to ê, Corollary 12 implies that y 0 and y 1 are paired in C to points in e that are located at distance at least 2r-ε from the extremities of e. As a consequence of Lemma 3 we have d Y (y 0 , y 1 ) ≤ (1 + 2ε/R)l(e) -4r + 3ε. Now, let (x, y), (x ′ , y ′ ) ∈ C such that x, x ′ ∈ e are two points at distance 11r/3 + 2ε from each endpoint of e. We deduce from Corollary 12 that y, y ′ are edge points, so d Y (y, y ′ ) ≤ d Y (y 0 , y 1 ) and from Lemma 3 that

d Y (y, y ′ ) ≥ d X (x, x ′ ) -(2d X (x, x ′ )/R + 1)ε ≥ l(e) 1 - 2ε R - 22r 3 -5ε
where for the last inequality we have used that l(e) = d X (x, x ′ ) + 22r/3 + 4ε ≥ d X (x, x ′ ). Putting all the above inequalities together we finally get

1 -κ(e) ≤ l(ê) l(e) ≤ 1 + κ ′ (e) with κ(e) = 10r 3l(e) + 5 l(e) + 2 R ε, κ ′ (e) = ( 3 l(e) + 2 R )ε
Using that l(e) ≥ b, we obtain that κ(e)

≤ κ = 10r 3b + ( 5 b + 2 R )ε and κ ′ (e) ≤ ( 3 b + 2 R )ε.
As a consequence, since e and ê are isometric to intervals, there exists a homeomorphism φ e : e → ê such that φ e is (1 + κ ′ )-Lipschitz and φ -1 is (1 -κ) -1 -Lipschitz. Since X and X are graphs, the homeomorphisms φ e can be glued all together to obtain a global homeomorphism φ : X → X such that φ is (1 + κ ′ )-Lipschitz and

φ -1 is (1 -κ) -1 -Lipschitz.
Recall that to each edge ê of X we assign a length equal to the diameter of the corresponding connected component in R 2r (E) plus 4r and we denote by d X the metric induced on X. To conclude the metric reconstruction part, we finally relate the metrics on Y and X. Theorem 14. There exists a map ψ : Y → X such that for any y, y ′ ∈ Y

(1 -κ) (1 - 2ε R )d Y (y, y ′ ) -ε ≤ d X(ψ(y), ψ(y ′ )) ≤ (1 + κ ′ ) (1 + 2ε R )d Y (y, y ′ ) + ε
with κ and κ ′ as in the Metric Reconstruction Theorem 6.

Proof. Let C be an (ε, R)-correspondence between Y and X. From the definition of correspondence, there exists a map (not necessarily continuous) f : Y → X such that for any y ∈ Y, (f (y), y) ∈ C. It immediately follows from Lemma 3 and Theorem 6 that ψ = φ • f verifies the desired inequalities.

Although the above result does not provide an explicit map, we provide an easy to compute map ψ : Y → X that satisfies similar inequalities when restricted to edge components. First we define ψ on the branch points: each branch point is mapped to the vertex of X corresponding to the connected component of R 2r (V) that contains it. We then define ψ on each connected component of R 2r (E). Let ê be an edge of X and let y 0 , y 1 be two points in the corresponding connected component in R 2r (E) such that d Y (y 0 , y 1 ) is equal to the diameter of this component. We parametrize isometrically ê by the interval [0, l(ê)]. Recall that l(ê) = d Y (y 0 , y 1 ) + 4r, we let ψ(y 0 ) = 2r and ψ(y 1 ) = l(ê) -2r. Now if y ∈ Y is in the same connected component of R 2r (E) as y 0 and y 1 we define ψ(y) = 2r + d Y (y, y 0 )

d Y (y0,y1) d Y (y,y0)+d Y (y,y1) . Lemma 15. For i = 0, 1 (1 -εM )d Y (y, y i ) ≤ ψ(y) -ψ(y i ) ≤ d Y (y, y i ) where M = 6/R + 1/b.
Proof. The proof of the case i = 0 and i = 1 being similar we give it for i = 0. Remark that ψ(y)-ψ(y 0 ) = d Y (y, y 0 ) (y,y1) and the second inequality is just the triangle inequality. Let x, x 0 , x 1 ∈ X be such that (x, y), (x 0 , y 0 ), (x 1 , y 1 ) ∈ C. Note that x, x 0 and x 1 are in the same edge of X so that d X (x 0 , x 1 ) can be expressed either as a sum or as a difference of d X (x 0 , x) and d X (x, x 1 ). Applying Lemma 3 three times and using that d Y (y 0 , y) ≤ d Y (y 0 , y 1 ) and d Y (y, y 1 ) ≤ d Y (y 0 , y 1 ) we obtain

d Y (y0,y1) d Y (y,y0)+d Y
d Y (y, y 0 ) + d Y (y, y 1 ) ≤ d Y (y 0 , y 1 ) 1 + 6ε R + ε Using that b ≤ d Y (y 0 , y 1 ) we finally get d Y (y,y0)+d Y (y,y1) d Y (y0,y1) ≤ 1 + ( 6 R + 1 b )ε.
From Lemma 15, we easily get the following corollary controlling the distortion on the metric induced by the restriction of ψ to the vertices of a connected component of R 2r (E).

Experiments

We implemented our algorithm in C++ using the Boost Graph Library [START_REF] Siek | Boost graph library[END_REF]. Experiments were conducted on a 2.33GHz Macbook Pro with 3GB of RAM. To assess the generality of our algorithm, we used four very different real world data sets: earthquake data, GPS traces, astronomical data and Image Webs. Table 1: Our algorithm was used on several data sets to reconstruct a simpler metric graph approximating the distances in the original graph. We randomly selected a sample of 100 points and computed all pairwise distances between points in the same connected components. The graph computation time is the total time of estimating degrees of nodes and reconstructing the graph. The original computation time shows the total time of computing these distances using the original graph. The approximate computation time is the total time it took to compute approximate distances with the help of the reconstructed graph. All times are in seconds.

Data Sets. We used four different data sets for which we expect there to be an underlying metric graph approximation. The first data set is that of earthquake locations through which we wish to learn topological and geometric information about earthquake faults. The raw data was obtained from USGS Earthquake Search [START_REF]Earthquake search[END_REF] and consists of earthquakes between 01/01/1970 and 01/01/2010, of magnitude greater than 5.0, and of location in the rectangular area between latitudes -75 degrees and 75 degrees and longitude between -170 degrees and 10 degrees. The underlying metric graph for this data set is the network of fault lines. The second data set is that of 500 GPS traces tagged "Moscow" from OpenStreetMap [2]. Since cars move on roads, we expect the locations of cars to provide information about the metric graph structure of the Moscow road network. The third data set consists of locations of galaxies in a portion of 3D space and there have been recent studies on the existence of filamentary structure in the distribution of galaxies [START_REF] Choi | Tracing the filamentary structure of the galaxy distribution at z ∼ 0.8[END_REF]. Lastly, we include an Image Web [START_REF] Heath | Image webs: Computing and exploiting connectivity in image collections[END_REF] data set which is a collection of images, with similar regions linked together to form a graph structure. Dense image collections are often acquired by mobile entities, and thus naturally contain long linear and circular parts, joined together at branch points.

Preprocessing and Parameter Selection. We performed some preprocessing to transform the raw data into a metric space (Y, d Y ) on which we could use our algorithm to discover a much simpler metric graph (X, d X ) approximating this space. Since real world data sets vary widely in both noise and scale, the specific preprocessing steps differ across the data sets. However, in most of our examples, we first construct a neighborhood graph on the data, and then used the shortest path metric space on the neighborhood graph as the input to our algorithm. The raw earthquake data set contains the coordinates of the epicenters of 12790 earthquakes in the latitude/longitude rectangle [-75, 75] × [-170 , 10]. As it contains outliers, we first preprocessed the data by removing points located in low density areas using the distance-to-measure function [START_REF] Chazal | Geometric Inference for Measures based on Distance Functions[END_REF]: points with average squared distance to their 30 nearest neighbors larger than 4.7 2 = 22.09 were discarded, resulting in the elimination of 284 points. Among the remaining data, points with average squared distance to their 50 nearest neighbors larger than 81 were also discarded to get a cleaner data set (eliminating 41 more points). Then, we randomly sampled 1600 landmarks among the points with average squared distance to their 50 nearest neighbors in the cleaned data set smaller than 1.5. Finally, we computed an α-complex with α = 4 on these landmarks, and used the shortest path metric on this complex as the input to our algorithm. For the road network data set, we first selected a metric ε-net on the raw GPS locations with ε = 5 using furthest point sampling. Then, we computed an α-complex on the ε-net as the neighborhood graph, but with α = 50. The astronomical data is similar to the earthquake data in the sense that it contains a lot of noise, which hides the filament structure. We built the input neighborhood graph on a set of landmarks selected in a similar fashion as that for the earthquake data set. The Image Web data set differs from the rest in that the raw data is a neighborhood graph, so no preprocessing was done. Our algorithm is parameterized by the spherical shell inner radius r, which in the analysis is allowed to be in a range of values that depends on a constant b that is the property of the underlying metric graph and the level of approximation attained by the data. In practice, however, we do not have an oracle for these constants. However, regardless of whether the assumptions in the analysis are satisfied, our implementation outputs a metric graph ( X, d X) and a map φ from the raw data to the metric graph. Using random sampling, we can estimate the level of metric distortion using X and φ. Thus, we are able to select the parameter by running our algorithm using various values of r, and checking for a balance between metric distortion and reduction of graph size. We note that even though the assumptions in the analysis may not be strictly satisfied, our algorithm returns a metric graph approximation that is in some cases dramatically smaller than the original data, while approximately preserving distances. In addition, we also varied the outer radius (5r/3 in the analysis) and the Rips-Vietoris parameter (4r/3 in the analysis) using the same process. Indeed, the constants 4r/3 and 5r/3 were chosen for ease of analysis. In particular, they ensured that in the proof of Theorem 8, all connected components were cliques, but in practice they may not be the best constants to use.

Implementation and Results

. Real world data sets often do not satisfy the assumptions we require for complete reconstruction, so we only replace connected components of R 2r (E) with edges of X if they are adjacent to exactly one (in the case of a self-loop) or two connected components of R 2r (V). Note that this process is local and hence it is possible to iterate this process in order to discover stratified structure at multiple scales. We also computed a map ψ from the original points to the reconstructed space X as described in Lemma 15. To evaluate the quality of the reconstructed graph for each data set, we randomly selected 100 points from the data set, and computed both original pairwise distances, and pairwise distances on X using ψ. We also evaluated the use of X to speed up distance computations by showing reductions in computation time. Statistics for the size of the reconstructed graph, error of approximate distances, and reduction in computation time are given in Table 1. Only pairs of vertices in the same connected component are included because we obtain zero error for the pairs of vertices that are not. We used these statistics to select the parameter r, as well as the outer radius of the spherical shell, and the Rips-Vietoris parameter.

The result of our algorithm on the earthquake data set is shown in Fig. 2. We observe two spurious branch points being detected on the component to the right as a result of the small stub sticking out between them. Nevertheless, our algorithm is able to replace the data by a much smaller graph, while maintaining small distortion of distance. Note that a trivial postprocessing step that removes all vertices of degree 2 could take care of the two spurious branch points. The GPS trace data set, shown in Fig. 3, provides the best results of all four data sets, showing a dramatic reduction in graph size along with a very small distortion of distance. This is expected considering that cars in most cities necessarily follow a road network, which fits the model of a metric graph very well. The metric graph structure in the astronomical data set, shown in Fig. 4, is much less apparent than that of the previous examples, and hence we were only able to reduce the graph size by one half. However, by doing so, we still approximately preserved distances and reduced distance computation time by more than 51%. The Image Web, shown in Fig. 5, was a very small example, and therefore suffers from metric distortion problems as noise levels are relatively large when compared to the size of the branching structures, but our algorithm was still able to reduce the already small graph size by 79% while keeping the median distance distortion below 18%.

Conclusion and Future Work

In this paper, we presented a first attempt at reconstructing a metric space of mixed dimension. We presented an algorithm with guarantees for the case of a metric graph, or equivalently, a 1-D stratified space. The same algorithm can be used to simplify the representation of a metric space that might possibly have an underlying metric graph structure. We also showed that, on real world data that doesn't perfectly satisfy the hypotheses, our algorithm still gives sensible and useful results.

A natural extension of this work would be to consider stratified spaces of higher dimension, as well as considering the data at multiple scales. Currently, we rely on the fact that our algorithm is relatively fast, and thus trying various scale parameters and checking for small reconstructed metric graph with small distortion is feasible. However, it is also interesting to consider the automatic selection of scales for which the data can be viewed as a reasonable approximation of a metric graph. We have also begun preliminary experiments for a multiscale version of our algorithm, which follows naturally from our implementation. It would be of interest to consider models of data where such a reconstruction gives theoretical guarantees. Other directions for further research include investigating the possibility of improving the distortion of the metric by allowing the addition of branch points to split edges that have too much distortion or to contract large regions of branch points into several points instead of just one. Having these options not only gives the user some choice on the trade-off between the size of the graph and the distortion, but also fits well with a multiscale approach.

Figure 1 :

 1 Figure 1: (a) A metric graph (in black) and 2 intrinsic balls (in blue and red). (b) Using a spherical shell to infer the degree of a vertex.

3 :

 3 Sy ← B Y (y, 5r/3) \ B Y (y, r) 4: degr(y) ← Number of connected components of Rips-Vietoris graph R 4r/3 (Sy) 5:

Definition 9 .

 9 Let (Y, d Y ) be an (ε, R)-approximation of X. Given 0 < r < R/2 and 0 ≤ α < 1, the (r, α)-degree deg r,α (y) of a point y ∈ Y is the number of connected components of the Rips-Vietoris graph with parameter 2αr and vertex set B Y (y, (1 + α)r) \ B Y (y, r) where B Y (y, r) denotes the intrinsic (closed) ball in Y with center y and radius r. Theorem 10 (Degree Inference Theorem).

Corollary 16 .

 16 If y, y ′ are in the same connected component of R 2r (E) corresponding to an edge ê in X then d Y (y, y ′ ) -εM l(ê) ≤ ψ(y) -ψ(y ′ ) ≤ d Y (y, y ′ ) + εM l(ê).

Figure 2 :

 2 Figure 2: Earthquake data: the input neighborhood graph is shown in cyan, the points marked as belonging to a branch are shown in red, and the points marked as belonging to an edge are shown in blue. The reconstructed graph is shown in dark blue.

Figure 3 :

 3 Figure 3: GPS traces: the input neighborhood graph is shown in cyan, the points marked as belonging to a branch are shown in red, and the points marked as belonging to an edge are shown in blue. The reconstructed graph is shown in dark blue.

Figure 4 :

 4 Figure 4: Astronomical data: the input neighborhood graph is shown in green and the reconstructed graph is shown in dark blue.

Figure 5 :

 5 Figure 5: Image Web: the input neighborhood graph is shown in cyan, the points marked as belonging to a branch are shown in red, and the points marked as belonging to an edge are shown in blue. The reconstructed graph is shown in dark blue.

Table 1

 1 summarizes our results, and a detailed discussion follows in this section.

		Earthquake	GPS Traces	Astronomical	Image Webs
	Number of Original Vertices	1600	28434	9276	530
	Number of Reconstructed Vertices	18	497	3651	112
	Number of Original Edges	3983	41669	34890	1711
	Number of Reconstructed Edges	9	5402	14808	409
	Graph Reconstruction Time	5.2846	43.2249	14.0829	0.729667
	Original Dist Comp Time	0.016386	0.777398	0.60322	0.021817
	Approx Dist Comp Time	0.004696	0.029821	0.29148	0.013379
	Dist Comp Time Speedup	249%	2507%	107%	63%
	Mean Distortion	6.4%	2.4%	22%	27%
	Median Distortion	8.8%	2.0%	19%	17%

see[START_REF] Gromov | Metric Structures for Riemannian and Non-Riemannian Spaces[END_REF] Chap.1 for the definition of the length of a continuous curve in a general metric space

We also include in our definition the 1-dimensional manifold isometric to a circle (one edge and no vertex)

Acknowledgements

This work was supported in part by NSF grants CCF 0634803, FODAVA 0808515, CCF 1011228, a grant from Google Inc., EU project CGLearning 255827, and ANR grant GIGA ANR-09-BLAN-0331-01. We would also like to thank the associated research team CoMeT.