An Integrated Ecological-Social Simulation Model of Farmer Decisions and Cropping System Performance in the Rolling Pampas (Argentina)

Changes in agricultural systems are a multi-causal process involving climate change, globalization and technological change. These complex interactions regulate the landscape transformation process by imposing landuseandcover change (LUCC)dynamics. Inorder tobetterunderstandand forecast theLUCCprocess we developed a spatially explicit agent-based model in the form of a Cellular Automata: the AgroDEVS model. The model was designed to project viable LUCC dynamics along with their associated economic and environmental changes. AgroDEVS is structured with behavioral rules and functions representing a) crop yields, b) weather conditions, c) economic profits, d) farmer preferences, e) adoption of technology levels and f) natural resource consumption based on embodied energy accounting. Using data from a typical location of the Pampa region (Argentina) for the period 1988-2015, simulation exercises showed that economic goals were achieved, on average, each 6 out of 10 years, but environmental thresholds were only achieved in 1.9 out of 10 years. In a set of 50-years simulations, LUCC patterns converge quickly towards the most profitable crop sequences, with no noticeable trade-o between economic and environmental conditions.

Introduction
. In agroecosystems, land-use and cover change (LUCC) is driven by simultaneous responses to economic opportunities, institutional factors and environmental constraints (Lambin et al. ). Di erent methodological approaches are used for assessing LUCC by capturing the dynamic process influenced by complex interactions between socio-economic drivers and biophysical conditions. For example, the spatially explicit LUCC models assess location suitability for di erent land uses and allocate changes to grid cells based on suitability maps (Verburg et al. a). Another modelling approach has been developed for capturing the behavior of the real actors of land-use change: individuals and/or institutions ("agents") become the objects of analysis, modelling and simulation, paying explicit attention to interactions between these agents of change (Castella & Verburg ). In this modelling approach, farmers' responses to environmental, economic and sociological constraints Region under study and cropping system description . AgroDEVS was applied to simulate LUCC dynamics in Pergamino, a typical location in the Rolling Pampas (Figure ), the most productive subregion of the Pampas where most of the annual cropping is concentrated (Hall et al. ). The Pampa is a fertile plain originally covered by grasslands, which during the s and s was transformed into an agricultural land mosaic by grazing and farming activities (Soriano et al. ). The predominant soils are typical Arguidols and the annual rainfall averaged mm (Moscatelli et al. ). AgroDEVS simulates LUCC using the most frequent crop types in the Pampa region (Manuel-Navarrete et al. ): ( ) the wheat/soybean double-cropping (W/S); ( ) maize cropping (M ), and spring soybean cropping (S). In this paper, the terms "crop"; "crop type" and "land use" (i.e. cropping systems) were used interchangeably as well as "farmer" and "agent". The agronomic decisions (i.e. genotype selection, fertilizer management, pest control, sowing date, and soil type) were representative of the most frequent situation for each of the cropping systems.

Cellular Automata-based modelling and simulation approach .
For the modelling activity we followed a scenario-based analysis, exploring agricultural mosaic dynamics. Then we adopted an agent-based modelling (ABM) approach, by identifying both landscape and agent-specific descriptors as parameters (fixed) or attributes (variable). We defined parameters as any fixed condition for describing the behavior or the condition of a model element. On the other hand, the attributes represent the system changes during the model run period. Lastly, we adopted a formal model-based simulation framework to specify mathematically both the parameters and attributes due to di erent behaviors. For this purpose, we encoded AgroDEVS using the Cell-DEVS formalism. Cell-DEVS is an extension for Cellular Automata of the more generic Discrete EVent System Specification (DEVS) formal modelling and simulation framework. On the one hand, the DEVS formalism permits to express and combine any kind of dynamical system (continuous, discrete event, discrete time) in a mathematical form that is independent of any programming language. On the other hand, Cell-DEVS provides the modeler with a meta-language tailored to facilitate the expressing of systems where the spatial arrangement of "cells", and their behavior, plays a key role. By adopting the DEVS-based approach, AgroDEVS becomes readily linkable with other DEVS models developed by others for di erent domains (e.g. climate, biology, sociology, and economics), potentially using heterogeneous techniques (e.g. di erential equations, equilibrium models, optimization models, stochastic processes). In Figure   .
Behavioral rules are used in the Cell-DEVS language to define the change of attributes that are local to each cell. These variables can express properties for that unit of space (e.g., in the -dimensional case it can be an agricultural plot). When the modeler uses such variables to express also the attributes of "agents" located at cells, then the ABM approach can be merged into the system. In the case of AgroDEVS each agent (namely, a farmer) remains fixed to his or her plot, and rules are used to express changes both for the physical environment and the for farmer, which eventually a ect each other. In this work, we adopted the CD++ simulation so ware toolkit, which is capable of interpreting Cell-DEVS and DEVS models and of simulating them. CD++ implements a generic DEVS "abstract simulator", consisting of a standardized algorithmic recipe that specifies how to simulate any DEVS model independently of any programming language of choice. We believe this generic, reusable and extensible DEVS-based approach provides AgroDEVS with very desirable scalability and sustainability features. The full model represents the collective (i.e., landscape-level) function that emerges from the aggregation of all farmers' outcomes. It also depends on exogenous variables (e.g., climate, crop (output) prices, and production costs) as well as endogenous variables (e.g., the farmer's technological level, the outcomes of neighboring farmers, and each farmer's performance history). The model proved particularly well suited to reproduce empirical situations where (a) there are changes in the relative production/output prices between potential land uses, (b) there is a specific climate regime that impacts on crop yields, or (c) there are varying aspiration levels for the farmers (Gotts & Polhill ). .
In Appendix E we provide more technical details about the model specification approach, its modelling language, rule specification, simulation framework and experimentation environment.

Model verification, validation and sensitivity analysis .
We submitted AgroDEVS to a systematic verification and validation process (Wilensky & Rand ). Firstly, a code walk-through (Stern ) was performed to review the model formulation and to ensure that all design concepts and specifications are correctly reflected in the code. In addition, AgroDEVS was run with very few farmers in the grid ( -) and results were examined closely (e.g., following dynamics of specific farmers, inspecting the outliers). The intrinsic complexities and uncertainties on both the magnitude and the nature of the forces drivingland use change, lead to expanding the scope of the straightforward evaluation between simulated and observed patterns in the model validation phase (Bert et al. ; Le et al. ; Nguyen & de Kok ). Our model development process seeks to validate the ABM through a process of continuous adaptation using feedback from the stakeholders (Ligtenberg et al.
). Specifically, AgroDEVS development entailed a continual discussion process with stakeholders (data not shown) from the study area in order to review both the rules and the assumptions of the model that initially came from the literature review. This approach directly engages stakeholders in model development by embedding it within the social process of policy development (Moss ). This changes the validation problem into an advantage: the agreement of participants or stakeholders may be an indicator of the validity of a simulation model (Troitzsch ). In AgroDEVS, the evaluation of the simulation is guided by the expectations, anticipations, and experience of the community that uses it for practical purposes (Ahrweiler & Gilbert ), and this supports the view that the very meaning of validity is dependent on the purpose of the simulation models under examination (Küppers & Lenhard ). Moreover, it is possible to develop a model that fits the data perfectly with a model structure that does not capture any real processes (Cooley & Solano ). Due to the complex nature of ABMs (e.g., nonlinear responses to parameters) a broad single-parameter SA (Appendix A) was performed (Railsback & Grimm ). This type of analysis assesses the e ect of each parameter over a wide range of values, as opposed to the traditional local analysis, in which each model input is varied by a standard small amount (Ligmann-Zielinska et al. ; Railsback & Grimm ). The goodness-of-fit of each simulated LUCC was calculated using ) Normalized RMSE: v = RM SE/Observed mean %T otal area; ) Probability of a Match: P M = #matches/(#matches + #mismatches); and ) Index of Observed Fit: IOF = (2 × P M ) − 1.
. P M are the chances that a simulation can correctly predict the order of any pair of observations. The number of matches (#matches) is calculated by counting the set of ordered pairs of observations that match the predicted ordered pairs of a simulation (Thorngate & Edmonds ). IOF is an index that derives from the calculations of the PM. Its values range from IOF = +1 (all observations match predictions), through IOF = 0 (half of the observations match), to IOF = −1 (none of the observations match predictions).

Model description
. The model description follows the ODD protocol "ODD" (Overview, Design concepts, and Details) to standardize the published descriptions of individual-based and agent-based models (Grimm et al. ). We describe in the following subsections ) the model's purpose, ) the system variables and ) the process overview to emphasize the main message of the model outputs. The remaining items of the ODD protocol (initialization conditions, submodels description, input data, modelling approach, and design concepts) can be found in the Appendixes B to F.

Model's purpose .
The main modelling purposes (see Edmonds ) of AgroDEVS are ) Prediction, as there is a need to be able to anticipate aspects of the agricultural system that are not currently known; and ) Explanation, as the model tries to identify the causal interactions between the main driving forces of LUCC phenomena (prediction alone does not provide insights on the internal dynamics driving the evolution of system processes).

System variables .
AgroDEVS maps farmers onto a regular grid in order to initialize the simulations. The model consists of two entities: ) the landscape and ) the agents that operate on the landscape. Each entity has its own set of fixed parameters and attributes that evolve throughout simulation cycles. The landscape parameters are a) the number of agents and b) the owner/tenant agent ratio and the landscape attribute is the overall outcome from the integration of all individual agent attributes results within the simulated landscape. The agent parameters are a) the land rental price (RP ) for the tenants, and b) its location on the grid. The attributes for each agent are a) the technological level (T L), b) the crop type allocation (or land use, LU ), c) the economic profit (P ), d) the renewability level (RL) of the embodied energy (emergy) consumption (see Appendix C Renewability level calculations section for an emergy concept explanation), e) the aspiration level (AL), f) the environmental threshold (ET ), and g) the weather growing condition (W GC).

Process overview .
AgroDEVS simulations advance with an annual time step, representing a single cropping cycle (CC) (Figure ).
Crop yields under di erent WGC were previously simulated (Appendix B, Table ) using the DSSAT model (Jones et al. ). We used local management data for defining resource level in each T L (e.g., fertilization regime, genotype) as well as five contrasting historical weather records for defining the di erent W GC levels. attributes and attributes calculated during each CC, respectively. Square boxes denote calculations performed at each stage (for more details of each step see Appendix C). Rectangles represent goal fulfillments, diamonds are conditional rules, and parallelograms are an agent's decisions. Solid and dotted lines denote Control flow and Data flow, respectively. Red borders highlight when the agent performs a neighborhood analysis, and grey background indicates that the calculated value shall be used in the next cropping cycle. .
At the start of the AgroDEVS simulation process (Figure , Stage ), each initial agent's configuration (LU t ) is exposed to a climate-related growing condition level W GC that ranges from very favorable to very unfavorable for high crop yield achievement. Then, based on ) the previously simulated crop yields, ) the crop price, ) the production cost of each T L, ) the land rental price and ) the crop type allocation into the farm household, AgroDEVS calculates the farmer's P t : the profit P for the current time t, or similarly, for the current cropping cycle CC. At this first step, AgroDEVS also calculates RL t (see Appendix C for RL calculation details) as a mea-sure of the renewable energy consumption of the cropping cycle. At the next stage, the initial agent's aspiration level (AL t ) is adjusted by means of the current W GC t level resulting in a new climate-adjusted aspiration level (CAL t ) (Figure , Stage ). Then, AgroDEVS uses P t and RL t values calculated in Stage for assessing the fulfillment of both the environmental (ET ) and economic (AL) goals (Figure , Stage ) of each agent. An agent will keep its crop type allocation once the P t value is greater than or equal to its CAL t at each CC, while nonfulfillment of environmental threshold (ET ) does not alter its crop type allocation decision. AgroDEVS assesses whether the farmer can upgrade or downgrade its T L for the next CC, according to their economic performance in the current CC (Figure , Stage ). For this adjustment, AgroDEVS defines a set of di erent working capital thresholds (W CT ) in order to access di erent T L values. The W CT fulfillment is assessed using P t . If P t is lower than the respective W CT , the farmer must lower its T L (and vice versa). Eventually, agents can remain at the lowest T L, regardless of its P t value (i.e., no agent is forced out of business). During this stage, AgroDEVS also adjusts the farmer's aspiration level (AL) for the next CC. This setting is based on a) the farmer's perception of the W GC level of the next CC, and b) the agent's failure or success at achieving the AL in the previous CC, respectively. Values of ) crop prices, ) adjustment factor of the aspiration level due to W GC, ) adjustment factor of the aspiration level due to T L, and ) working capital threshold for T L are shown in the Appendix B, Tables to , respectively).
The simulated LUCC patterns replicated the overall trend towards soybean-dominated landscapes observed in the region since the mid-s. The ordinal pattern analyses (OPA) showed a similar, relatively high P M for maize, wheat/soybean, and soybean, evidencing the model's capability to predict ordinal (higher/lower) values of crop type cover. However, the accuracy in predicting the magnitude of these changes was lower (Figure ). Based on the RMSE method, the goodness-of-fit of the wheat/soybean simulated cover was the lowest among the three-crop types, resulting in a significant overestimation of this cropping area. The model underestimated simulated soybean cover. In addition, the v values (i.e., the RMSE relative to the observed mean) was remarkably higher for wheat/soybean than in both maize and soybean. Profit, renewability, and technological levels . Profit (P ) and renewability level (RL) values for the Pergamino simulation correlated without any significant trend throughout the study period ( Figure ). Remarkably, in only four years of the simulated period, the simulated landscape exhibited higher RL values than the ecological threshold (ET > 50%RL). Inter-agent variability in the simulated landscape depicted a range of decisions (in terms of crop type allocation) that resulted in a great variability on average P values during the simulation ( Figure   . The interquartile P range ( -% percentile) varied between and US$/ha (Figure , Panel A). However, the RL interquartile range was not as high as the observed in P (Figure , Panel A). The relative homogeneity between each cropping production system, in terms of external inputs, generated a remarkably low variability in RL terms during the whole simulated period. In this case, the interquartile RL range varied from . % to . % with a median value of . %, almost ten percentage points lower than the environmental threshold (ET ) of % (Figure , Panel A). While in terms of RL, the variability between agents was significantly lower than those observed for P , it was possible to detect both maximum and minimum RL values showing the model capability to simulate a wide area of decisions during the simulated period. Minimum average value for an agent outcome in terms of RL was % and the maximum RL average value for an agent was . % (Figure , Panel A), a condition of very close agreement with the environmental threshold (ET ). Agent behavior in AgroDEVS was also assessed by inspecting the intra-agent variability, using both P and RL during the simulation period ( Figure , Panel B). The interannual coe icient of variation of P showed a median value of %, with an interquartile range between . % and . % ( Figure , Panel B). The RL interannual variation was significantly lower than P , exhibiting a median coe icient of variation of . % and an interquartile range between % and % ( Figure , Panel B).

Figure :
Inter-agent and intra-agent variability of the Pergamino -simulation, and the Ecological threshold (ET ) and economic aspiration level (AL) goal agreements. Panel A shows the average profit (P ) and renewability level (RL) of the agents in the Pergamino -simulation. The horizontal solid line shows the overall average of P and RL (n = 625). The extremes of the whiskers represent the % and % quartiles, and the numbers show the minimum and maximum P and RL average values for all agents. Panel B shows the coe icient of variation of Profit (cv P ) and of Renewability Level (cv RL) for all agents. The horizontal solid line shows the overall average of cv P and cv RLV (n = 625), the extremes of the whiskers represent the % and % quartiles, and the numbers denote the minimum and maximum cv P and cv RL values. Panel C shows the average percentage of agreement (% goal agreement) for all agents of the Pergamino simulation. The horizontal solid line shows the overall average of ET and AL goal agreement (n = 625), the extremes of the whiskers represent the % and % quartiles, and the numbers show the minimum and maximum ET and AL goal agreement values. .
An assessment of the simulation results was carried for both the environmental (ET ) and economic (AL) goal agreements exhibited by agents throughout the simulation period ( Figure , Panel C). The goal agreement metric indicates the percentage of years during which an agent fulfills each of the goals (environmental and economic). The agents exhibited an interquartile range of economic (AL) goal agreement between % and % with a median of % and maximum and minimum of . % and . %, respectively (Figure , Panel C). These values are noticeably higher than the RL goal agreement that showed a median of . % with an interquartile range between % and % (Figure , Panel C). The agent's capability for adjusting the AL, based on both the current climate condition (W GC) and the AL fulfillment in the previous agricultural cycle, could explain the better goal agreement when compared against RL, which is a fixed parameter. When inspecting the extreme values for RL goal agreement, it was possible to identify a single agent exhibiting a maximum value of . %. This means that, under the same economic and climate conditions, this agent was able to fulfill nearly out of years the RL threshold (ET ) by means of its crop type allocation decisions.
. Lastly, the model outcome was also assessed by inspecting the final distribution of the three di erent technological levels (T L) across all agents of the Pergamino simulation ( Figure ). The closely related observed and simulated agent distribution underline the model's ability to represent the process of agricultural intensification evidenced by a higher proportion of farmers using a technological level with higher levels of inputs. Both observed and simulated T L distribution patterns exhibit a dominance of the high T L. However, the model retains a greater percentage of agents at the lowest level (T L low) than in the observed data. It is possible that the model structure, which allows agents to remain using the lowest T L despite not reaching the minimum wealth to face those costs, is an explaining factor for the overestimation of Low T L at the end of the simulation. The objective of the long-term simulations was to assess the magnitude of the e ect of the tenure condition (i.e., an owner or tenant-dominated landscape) and climate (i.e., five contrasting climate regimes) on both LUCC (  Profit and Renewability level . The variability between long-term scenarios was lower for both P and RL than for the simulated LUCC ( Figure  ). Regarding the variability induced by the climate regime, the simulated landscapes were stabilized at increasing P -values, as climate scenarios were better (Figure : Panels L, A, and H). This e ect of profit improvement occurred in both tenant-and owner-dominated landscapes. In this latter case, the stabilized P values were due to the di erential income associated with non-payment of land rental. .
Although the simulated landscape configurations were clearly di erent ( Figure ), the RL variability among long-term scenarios under constant climate regimes (when stabilized) showed very small changes (less than %) between the maximum and minimum RL final values (

Discussion
. Determining patterns of land use change is undoubtedly a multifaceted challenge. When relying on simulations, several factors come into play such as climate influence on crop yield, the farmers' decisions, economic prices and costs, and the cognitive description of farmer behaviors (Hare & Deadman ). Thus, the construction of a LUCC simulation model entails necessarily the coupling of social and environmental models (Müller et al. ). In this paper, we developed an ABM expressing several of these characteristics. AgroDEVS simulations were able to reproduce observable LUCC trends of the most representative cropping systems in the region under study. An ABM validation process has the peculiarity of being subjected to conflicts between achieving accuracy in matching the outcome of a simulation or in the processes simulated (Brown et al. ). This trade-o is usually solved based on the research goals. In the case of LUCC simulation models, both aspects are important. Predicting the LUCC trends is extremely important in decision making by policy makers (Verburg et al. b). However, this information should be supplemented through an understanding of the underlying LUCC processes involved. This is required to identify potentially unsustainable land use regimes and correct them. AgroDEVS' results show that its structure is able to detect the overall trend on land use changes. This is done through clear and explicit modelling that reflects key process dynamics such as the climate influence on crop yields, farmer decisions and the landscape emergent properties due to farmer interactions at smaller scales. .
The contribution of the model developed in this work can be assessed by analyzing carefully the LUCC simulation results for Pergamino -, while considering the trade-o s between output accuracy and processes understanding. Although the model was able to simulate the land use change dynamics of the three crops analyzed (i.e., the ordinal fit is always higher than . ), the adjustment based on the distance between observed and predicted (i.e., the v value) could be improved by the inclusion of other variables or exploratory processes (e.g., agricultural policy decisions not considered; the dynamics of prices and costs, etc.). However, the current model structure maintains the relative profit between activities which is highly sensitive to environmental conditions (W GC). Thus, the distance between the observed and predicted LUCC can be used as a predictor of the di erence between what could have been a LUCC trend (based solely on the response to the environment in order to maximize profit, i.e., LUCC simulated) and another path that did not strictly follow the parameters of maximizing profit (i.e., LUCC observed). Notably, the region studied has been frequently subjected to decisions in agricultural policy (e.g., di icult marketing of some crops, imposition of export taxes) that strongly influence the changes in land use through mechanisms not directly related to profit activity (Porto & Lodola ).
. From the formal modelling and simulation point of view, the DEVS approach (and its related Cell-DEVS spatially specific flavor) o ered several salient features, already discussed in the model description section. Notably, in the context of socio-environmental systems, there is a feature that stands out. Namely, the input-output portbased hierarchical modularity permits the design of interdisciplinary models by composing complex systems through the interconnections of simpler ones. For instance, in AgroDEVS, the DEVS Atomic Model representing climate dynamics can be replaced by another, more accurate or sophisticated DEVS model developed by experts in the climate domain, without requiring to alter the landscape (Cell-DEVS) portion of the system. This approach is consistent with a current trend towards Systems of Systems-oriented modelling and simulation (Zeigler et al. ). This is especially relevant for socio-natural-economic systems, whose domain-specific submodels are constantly subject to revisions, improvements or replacements. Regarding the use of Cell-DEVS to spatially represent farms, abstract cellular landscapes have already been used in other ABM-LUCC models (e.g. FEARLUS, AgroPoliS, Pampas Model; Polhill et al. ; Happe et al. ; Bert et al. ). Results from AgroDEVS indicate that this approach has proven to be an e ective and simple way of spatially modeling agricultural systems. .
Incorporating environmental assessment in LUCC simulation models is a very desirable feature that is beginning to be explored (Veldkamp & Verburg ). The analysis of environmental impacts on managed ecosystems has o en been applied based reductionist approaches, identifying changes very accurately, but reducing the relevance of the results by not addressing an integrated or holistic approach for ecosystems modelling (Shanmuganathan et al.
). The AgroDEVS structure acknowledges the need for a systemic modelling by including an emergy (embodied energy) renewability level as a proxy for assessing ecosystem functioning of the cropping systems studied, an approach that has been tested before with agricultural systems both in the studied region and in other agricultural ecosystems (Agostinho et al. ; Dewulf et al. ; Ferraro & Benzi ). Concerning the use of environmental work, the simulated ecosystems did not show a clear dynamic of increase or decrease in its reliance on energy from the economic system. Moreover, the results from the simulations, both in the Pergamino -and the long-term runs, also failed to show a clear trade-o between the environmental performance (assessed through renewability level) and economic performance (assessed through economic profit) in the studied systems. The simulated dynamics showed that the systems reduced their dependency on exogenous emergy (i.e., fossil energy) to the extent that environmental conditions improved, and thus natural resources become more important. Clearly, this happens because at increasing levels of technology, increased use of external emergy is roughly proportional to the increase in capturing endogenous emergy associated with growth conditions (W GC) that improve crop yields. The high emergy return seems a characteristic of the farming systems of the Pampas region at the field scale (Ferraro & Benzi , ); but the AgroDEVS simulations showed that this pattern is maintained at the landscape scale due to the emergent properties that arise from the integration of individual farmers behavior. .
In the context of changes in agricultural policy in the area under study, the results of AgroDEVS seem to indicate the need for policy options that alter the relative prices of crops, to avoid the predominance of monocultures or systems highly dependent on the soybean crop. In this context, the AgroDEVS capability for identifying potential trade-o s between di erent agroecosystem domains (i.e., economic, environmental) is extremely important in the diagnosis of agricultural sustainability (Tittonell ). Policy options, which have direct e ects on model variables such as prices, can be readily tested into AgroDEVS for their repercussions on farm incomes and LUCC emergent patterns. Moreover, there is always the chance to include other variables in the modelling exercise. In the interest of better representing the heterogeneity amongst di erent farmer's decision logics in AgroDEVS, future e orts could be aimed at exploring the farmers' decision-making process. This exercise could expose new relevant variables that improve the representation of agents' behavior. However, this inclusion would require a new setting and a more complex numerical validation. The cost-benefit balance of these additions should be analyzed carefully to avoid incurring in an overfitting, conspiring against the understanding of the true underlying phenomena under study (Brown et al. ).

Conclusion
. The model presented in this work represents an e ort to integrate, within a complex system simulation, the e ects of weather on crops, the farmer decision logic, and the cropping system profit as main LUCC drivers. AgroDEVS simulations based on real landscape data showed that LUCC direction was better represented than its magnitude, in terms of the land cover dynamics. The long-term simulations showed a dominance of cropping systems that included soybean crop, and this dominance was stronger for monospecific soybean crop in scenarios under constant climate. The double cropping W/S dominated mainly in scenarios under a variable climate. When assessing the farmer condition e ect (i.e., tenure and climate) on LUCC, tenure resulted in much less e ect on LUCC than the weather conditions (W GC). Finally, simulations showed no trade-o s between environmental and economic outcome both in simulation used to validate the model and long-term scenarios. The results suggest that LUCC modelling and its environmental and economic consequences is feasible and useful using an ABM approach. The AgroDEVS simulations allow not only to speculate on the LUCC dynamics, but to gain a greater understanding of the underlying processes involved. Future research should focus on improving the model structure to include di erent agent behaviors (e.g., multiple agent's profiles) as well as social and political factors both for predicting the LUCC direction and to assess their relative magnitudes more accurately. Moreover, incorporating di erent productive regions of Argentina into AgroDEVS could expand the conclusions achieved in this paper, revealing di erences in LUCC dynamics at an ecorregional level.
ulator is CD++, a general purpose discrete-event simulator written in C++ by several generations of developers. The version recommended for this work is: https://github.com/SimulationEverywhere/CDPP_ExtendedState s-codename-Santi.To access the model please contact the authors.

Appendix A. Sensitivity Analysis
In a long-term simulation ( years), we tested the e ect of ) crop prices, ) W GC, ) the Owner/Tenant agent ratio and ) rental price. We run AgroDEVS varying one parameter at a time over a range of values (Table ), while keeping the other parameters at their reference value (i.e. the initialization conditions of the Pergamino simulation; for more details see Appendix B). The range of values for the crop prices was extracted from historical data and are the lowest and highest price of each crop for the years -.

Sensitivity analysis
The WGC scenarios were created by alternating % of the campaigns with historical values (i.e. actual weather for that year) and the other % with the WGC level desired. Thus, five scenarios were created comprising % of historical values and % W GC (V U , U , R, F and V F ). The owner/tenant agent relation ranged from % of the agents being owners and % tenants to % owners and % tenants in increments of % (i.e. / , / , / , / , / ). The rental price ranged from US$/ha , to US$/ha , (Table ).
The sensitivity analysis showed that the model is sensitive to wheat and soybean price changes in terms of profit. When the price of one of these crops rises, the land use cover that dominates the grid is W/S. This crop sequence has a high profit potential, and thus explains the change in profit observed. However, if the soybean price drops below the reference value the farmers cannot compensate with another LU and the profit drops. In terms of renewability, the model was not sensitive to changes in crop prices and all crops behaved similarly. Only a large increase in wheat price caused a % reduction in the RL (Figure ).  The model's response to di erent weather patterns was clear in terms of change in profit (Figure , le ). Under favorable W GC, the profit was higher than the reference value. In the best-case scenario ( % of the years under V F weather), the profit was increased by % of the reference value. Conversely, under unfavorable W GC the profit tended to drop. In the worst-case scenario, ( % of the years under V U weather) the profit decreased %. This behavior is explained by the quick dominance of a more profitable LU (W/S) and management level (High) under favorable weather conditions, and the dominance of a defensive LU such as soybean under unfavorable weather conditions. The model showed little response in terms of RL under di erent W GC (Figure ,  right). The results indicate that under favorable weather conditions the RL could drop slightly due to the dominance of W/S. On the other hand, unfavorable weather conditions would lead to an increase in S, provoking an increase in RL. The model was slightly sensitive to changes in rental cost, leading to profit variation (Figure , le ). As expected, an increase in rental cost diminished the profit in tenants resulting in a decrease in the average landscape profit (no dominance of LU was observed here). On the other hand, the RL was not sensitive to changes in rental cost (Figure , right). This could be explained by the relatively low fraction of the profit that the rental cost represents.

Appendix B. AgroDEVS Initialization Conditions
Crop yields are based on simulations using crop models in the Decision Support System for Agrotechnology Transfer (DSSAT) package (Jones et al. ) that have been calibrated for the studied location (Mercau et al. ). Crop management of each T L (e.g. genotype, sowing date, fertilizer rate, pesticide use, tillage operations, and production costs) was used for running the DDSAT crop yield simulations, and was based on local management data. Crop yield variability during the DSSAT simulated period ( -) was used to obtain the crop yield under contrasting W GC levels. The five W GC levels (very unfavorable, unfavorable, average, favorable, and very favorable) correspond to di erent percentiles (i.e. , , , , and , respectively) of yields simulated using historical weather records. DSSAT simulations do not account for important factors such JASSS, ( ) , http://jasss.soc.surrey.ac.uk/ / / .html Doi: . /jasss. as weeds, diseases, and pests. Thus, we empirically adjusted the attainable crop yield (van Ittersum & Rabbinge ) resulting from DSSAT simulations in order to model actual crop yield. Local data of simulated versus observed crop yield were used for obtaining the adjusting coe icients (attainable to actual yield) at each T L for each crop species (Mercau et al. , ; Satorre et al. ). AgroDEVS reflects three incremental T L: low (L), average (A) and high (H) technological levels; and three crop types (maize, soybean and wheat/soybean double-cropping, see region under study and cropping system description section for production system description). Using DSSAT and historical records of production costs, two of the look-up tables are built: ) a crop yield matrix (Table ) and ) a cost matrix for representing the full crop type ( levels), T L ( levels) and W GC ( levels) combination (Table ). As for crop yields and production costs, the emergy accounting method is used for building a third look-up table for representing the full combination of crop types ( levels), T L ( levels) and W GC ( levels) ( Table ). Historical prices of maize, soybeans, and wheat, as well as production costs (i.e. fertilizers, seeds, pesticides, and harvest and sale costs), were extracted from the Argentine trade magazine "Márgenes Agropecuarios" (http://www.margenes.com). In all scenarios (i.e. the Pergamino simulation and the long-term simulations) we assumed constant output prices equal to median for -.

Maize Soybean
Wheat/ Soybean Table : Adjustment factor αAL(W GC) of the aspiration level (AL) due to weather growing condition (W GC) level. W GC classes: Very Unfavorable; Unfavorable; Average; Favorable; Very Favorable.

Very Unfavorable
Unfavorable Average Favorable Very Favorable -.
-. . . . Table : Adjustment factor αAL(BN ) of the aspiration level (AL) due to agent technological level (T L) and the best neighbor technological level T L t (BN). A positive (or negative) αAL(BN ) value indicates an increase (or decrease) in AL since the best neighbor exhibits a higher (or lower) T L than that evaluated by the agent. -.
-. . Table : Working capital threshold (W CT ) for each technological level (T L) class (expressed in US$/ha.). W CT values equal to % of the highest production cost for each T L, considering and average indebtedness rate of % (AACREA ). T L classes: H (High); (A) Average; L (Low).

Appendix C. Submodels Description
Profit calculations P is calculated as the gross income (yield times product price) minus direct costs. Direct production costs include fixed and variable components. Fixed direct costs do not depend on crop yield (e.g., seed and agrochemicals). Oppositely, variable direct costs are yield-dependent (e.g., harvest, marketing fees and grain transportation). Finally, each farmer P is calculated using the individual crop profit a ected by each crop type allocation into the agent farm. Figure is an excerpt of the rule that calculates P : an agent has a defined LU and T L and is capable of sensing the WGC, and therefore it can calculate its P by multiplying its crop type allocation, yields and prices, and subtracting costs. Figure : Excerpt of the rule that performs agents' profit calculations. In the cases of yield, price and cost, #macro involves lookup tables.

Renewability level calculations
AgroDEVS calculates Renewability Level (RL) values by using the emergy synthesis procedure (Odum ). Briefly, the emergy accounting methodology tries to account for both the natural and human-made capital storages. The emergy accounting method values these storages using a common unit of reference, the solar equivalent joule (seJ). The method accounts for the environmental support provided directly and indirectly by nature to resource generation and processing; it focuses on valuation of the intrinsic properties of ecosystems (Mellino et al. ). For further details on emergy synthesis methodology see Brown et al. ( ) and Ferraro & Benzi ( ). AgroDEVS uses the renewability level as a sustainability metric (Giannetti et al. ). RL is then calculated as the ratio of renewable emergy to total emergy use, as it follows: where, RL(%): renewability level value, R(seJ/ha −1 year −1 ): renewable flows from nature, N (seJ/ha −1 year −1 ): nonrenewable flows from nature, F (seJ/ha −1 year −1 ): imported economic flows, S(seJ/ha −1 year −1 ): services.
Therefore, at each model time step (i.e. a CC) each individual farmer renewability level (RL) is used as a measure of environmental impact. In the figure below ( Figure ), an agent uses its LU and T L along with the sensed WGC to obtain the final farmer's RL value (by a ecting the individual crop renewability level with the crop type allocation into the agent's farm).

Figure :
Excerpt of the rule that performs agents' renewability calculations. In this case, the #macro(emergy) involves emergy lookup tables.

Agent decisions
Both P and RL values for each agent are related to the AL and ET thresholds, in order to trigger the agent's decisions. In AgroDEVS, the economic goal (AL) is dynamic and it is based on the aspiration level adjustment. It represents the currently dominating economic paradigm, where the receiver (the agent) is the market actor who decides the system outcome value (Grönlund et al. ). Oppositely, the environmental threshold (ET ) is fixed and it represents a strong sustainability view, where the value approach is grounded in systems science rather than economic science, where a value focused on the system level is accepted (Grönlund et al. ). During the model simulation process, the fulfillment of the economic goal drives the crop type allocation. Thus, the LUCC process is triggered when the economic threshold (i.e. aspiration level) is not accomplished.

Aspiration level adjustment
A first AL adjustment is based on the W GC, as it follows: where CAL t = climate-adjusted aspiration level (US$/ha), AL t = aspiration level (US$/ha) calculated at the end of the previous CC, αAL(W GC) = adjustment factor of AL due to W GC level for the current CC (see Appendix B, Table for αAL(W GC) values used in the simulations).
A second AL adjustment defines the aspiration level for the next CC and it is based on the learning and adaptation model (Bert et al. ). If P t > CAL t then the next aspiration level is subjected to an incremental adjustment using a weighted average, as it follows: where AL t+1 = aspiration level (US$/ha) for the next CC, CAL t = aspiration level (US$/ha) a er current CC adjustment, P t = Profit (US$/ha) calculated in profit calculations section. Figure : Excerpt of an agent's aspiration level adjustment. In this example, the agents' profit was larger than its AL, so it proceeds to perform an incremental adjustment.
In the figure above we can see an excerpt of the rule for adjusting incrementally the aspiration level using a weighted average when an agent's P t > CAL t . When both the agent and its neighbors finish with the calculation of their P , RL and the fulfillment of both economic and environmental goals, the agent proceeds to adjust its AL for the next cropping cycle. When the economic outcome is lower than the aspiration level (i.e. P t < CAL t ), the agent's perception is extended to include the influence of the physical (Moore) neighbors, and the AL t+1 is adopted by inspecting the neighbors' P outcomes. Thus, the farmer adopts AL t+1 using the neighbors' P data. If there is at least one neighbor under the condition of P t > CAL t then the farmer selects its best neighbor (BN ), in profit terms, and the AL t+1 is calculated using both the AL for the best neighbor, CAL t (BN ), as well as an adjustment factor due to di erences in T L between farmers: where: AL t+1 = aspiration level (US$/ha) for the next CC, BN = the agent with the highest P value in the Moore neighborhood (n = 8), CAL t (BN ) = aspiration level (US$/ha) for the BN in the current CC, αAL(BN ) = adjustment factor of AL due to BN technological level (see Appendix B, Table for αAL(W GC) values used in the simulations).
In the case that P BN ≥ CAL t (i.e. no neighbor meets the economic threshold) then the next aspiration level is subjected to a detrimental adjustment using a weighted average, as it follows: where: AL t+1 = aspiration level (US$/ha) for the next CC, CAL t = aspiration level (US$/ha) a er current CC adjustment, P t = Profit (US$/ha) calculated in Profit calculations.
Di erent values for incremental and detrimental adjustment ( . and . respectively) were applied in order to simulate the farmer willing to tolerate higher payo s more rapidly than lower ones, thus showing greater resistance to downward changes (Gilboa & Schmeidler ).

Technological level adjustment
Farmers are also able to upgrade or downgrade their technological level (T L) a er inspecting their own P outcome in the CC. The rules for adjusting the T L are based on the working capital threshold (W CT ) which represents the highest production cost within each T L ( Table ??). Rules for T L adjustment are: where, T L t+1 = Technological level for the next CC (H: high; A: average and L: low), W CT (T L i ): working capital threshold (average or high) for T L i (low, average and high).
In Figure we can see an excerpt of the rule for adjusting an agent's T L. When both the agent and its set of neighbors are done with calculating their P , RL and evaluating the fulfillment of both economic and environmental goals, the agent proceeds to adjust its T L for the next cropping cycle by comparing its P with the W CT . Figure : Excerpt of an agent's technological level adjustment. In this case, #macro (wc max ) involves WCT lookup table (Appendix B, Table ).

Land use and cover change (LUCC)
Final farmer decision defines the land use configuration. The LUCC process is triggered when the economic threshold is below the economic outcome and also there is at least one neighbor in the Moore neighborhood that meets the economic threshold as follows: where, P t = Profit (US$/ha) calculated as in profit calculations section, P t (BN ): Best Neighbor Profit (US$/ha), CAL t = Aspiration level (US$/ha) a er current CC adjustment.
A er inspecting this condition, the agent selects the crop type allocation of the best neighbor (BN ) as follows:  In the figure above there is an excerpt of the rule that copies the best neighbor's land use allocation. If both an agent and its neighbors have already calculated their P and RL then the model will proceed to select the JASSS, ( ) , http://jasss.soc.surrey.ac.uk/ / / .html Doi: . /jasss. neighbor with the highest P and copy its LU allocation, AL and T L to be used later in the simulation and, if necessary, assigned to the agent.

Pergamino -simulation
Input data from -were used for both the model initialization and to test the model outcomes against the actual land use changes in the studied area. Table presents a summary of these initial conditions.
The observed LUCC dynamics were obtained from agricultural surveys (SAGPyA ), and both the initial ( ) and final ( ) T L frequency distribution among farmers were obtained from the national agricultural census (BOLCER ; INDEC ). In order to represent the main tenure regimes in the studied area (owners and tenants), AgroDEVS was initialized using the Owner/Tenant relation (INDEC ). The land rental price was set to historical values of . t of soybean per hectare (Margenes Agropecuarios ). Crop yield and renewability level values under di erent weather conditions (W GC) and technological levels (T L) were calculated as it is explained in Appendix C profit calculations section. In the initial landscape configuration, each farmer was assigned an initial working capital (W C) according to his initial T L. Initial AL level were fixed for each farmer in order to account for the % of the direct costs (i.e. 0.6W CT ) of the T L adopted in the first CC of the simulation (AACREA ). Oppositely, the ET value remains constant during the simulations and is fixed in RL = 50% (Table ). The accuracy of the simulations was assessed using both ) the squared distances between a simulation's outputs and a set of observations (RMSE) and ) the ordinal pattern analyses (OPA) (Thorngate & Edmonds ). OPA indicates the topological fit between observed and simulated outputs but does not consider their closeness.

Long-term scenarios
AgroDEVS was also run over a -year period, under contrasting scenarios due both ) five W GC regimes (constant unfavorable, constant average, constant favorable, a see-saw pattern of very unfavorable-average-very favorable, and a random regime), and ) two tenure regimes based on the landscape pattern of Owner/Tenant agent relationship ( / and / ). As tenants are more focused on short-term income and are less likely to invest in longer-term management strategies than owners (Soule et al. ), we used AgroDEVS for testing the hypothesis that tenant farmers are less likely than owner-operators to adopt crop allocation decisions that lead to sustainable LUCC trajectories. All crop types and T L were set to equal distribution in the landscape (i.e. % of the total area for each crop and each T L) at the initialization. However, the internal assignment of each crop type for each agent was set randomly. Model simulations were inspected in terms of the dynamics of ) crop type allocation of total area, ) profit and ) renewability level.

Appendix E. Modelling Approach
In Figures and we provide an illustrative, non-exhaustive sample view of the Cell-DEVS specification used for AgroDEVS. Most of the information declaring the structure, components and behavior is found in a text file with .ma extension. Figure : Sample excerpt of an AgroDEVS.ma definition file (Le ) Schematic of the AgroDEVS system. The "climate" DEVS model is connected with the "landscape" Cell-DEVS model (Right).
In Figure we can see the main statements defining structure, components and interconnections through input-output ports. In Figure we provide an excerpt of the agrodevs.ma file highlighting relevant lines. Note that three dots ". . . " in the figures denote lines omitted for the sake of brevity. In Figure  we can see main statements defining behavior for the landscape cellular automaton. Within each cell, a list of rules is evaluated sequentially in a top-down fashion, using a Value, Delay, Condition structure: the first rule that evaluates its Condition to true, will apply the Value to its attributes a er a Delay amount of time. Once a rule is applied, the simulator recommences the cycle evaluating the rules from the top. This mechanism is applied asynchronously, simultaneously and in parallel to all cells in the model. The global timing for the whole cellular space is an emergent property driven by the local timing applied by each independent cell. The pairs (X, Y ) denote the positions for neighbors relative to each currently evaluated cell denoted with (0, 0).

Simulation framework and so ware experimentation environment
As discussed previously in the modelling approach section, AgroDEVS uses a DEVS-based formal approach. One salient feature of the DEVS formalism is the strict and clear separation between simulation algorithms and model specification. As a specialization of DEVS for Cellular Automata, the Cell-DEVS formalism inherits this separation. Di erent Cell-DEVS-capable simulators should be able to simulate a given Cell-DEVS model. In turn, di erent interactive user interfaces can be used to help with the design and maintenance of model specifications and to interact with the simulation exercises. Figure   The CD++ generic simulation engine for DEVS and Cell-DEVS models resides in a Simulation Server. It can retrieve and simulate Cell-DEVS formal specifications of models that are stored in a Database Server. Such models can be written directly by specialists in the Cell-DEVS formalism or by higher-level experts in the agricultural domain using a friendly web-based interface. Once models are defined, they can be reused to launch as many simulation experiments as required, by invoking the simulator also from a web-based experiment interface. Figure describes sequences of steps for two typical use cases: Model Specification and Simulation, denoting the components involved in each of them.

Objectives
Agents pursue to achieve a P level in each CC above their AL while satisfying (or not) a fixed RL. If P is below agents' AL, they will be unsatisfied and will seek a di erent crop type allocation in the agents' neighborhood. If an agent's capital availability drops below the W CT for each T L, they must adopt a lower-cost T L, but agents never quit farming. Both landscape and individual RL emerge from the crop allocation decided for each agent, but it is not used as a farmer s goal. Rather, this is an emergent property due to local rules.

Prediction
Agents predict the future consequences of their decisions (i.e., they build expectations) about P based on past outcomes and current weather information (W GC), following the theory of adaptive expectations (Shell & Stiglitz ). Agents revise their AL in each CC using adaptive rules based on a) the dissimilarity -in the previous CC-between the predicted and the observed P , and b) the seasonal climate forecasts (W GC) for the current CC (see sub-model section for algorithms adjusting AL).

Sensing
Agents have information about their capital availability and consider this variable in their decisions about the potential adoption of higher (or lower) T L. In addition, they are aware of the P achieved by their eight neighbors (Moore neighborhood) during each CC. Finally, farmers are informed about the expected status of external contextual factors (i.e. seasonal forecast) for the current CC.

Agent-Agent interaction
Agent interaction is based on the imitation of the crop type allocation of the best neighboring agent (i.e., the highest P achieved in the neighborhood) in the case of dissatisfaction with its own P achieved (as compared to its AL). Otherwise, the agent repeats its current crop type allocation in the next CC.

Agent-environment interaction
Crop yield simulations reflect the interaction between farm decisions and both climate and soil conditions. These models capture the T L e ect (through increasing fertilization and pesticide use) as well as the W GC of each CC (mainly global radiation and precipitation) and land quality factors (which are constant during the simulation). As crop simulation models represent biotic-unconstrained yield potential, final yields are adjusted in each combination of crop type-T L using potential weeds and pest infestation losses (see simulations input section for an example of crop yield data). Agent decisions have no e ect on resource dynamics, through crop type allocation, and it is only on the current W GC that agents receive feedback about changes in resource conditions.

Stochasticity
Stochasticity is used to assign the T L and the crop type allocation to each agent only at the initialization procedure.

Collectives
Agents do not form or belong to aggregations that a ect or are a ected by other agents in the model version described here. Further model refinements should include collectives.


Testing for error invariance in separable instrumental variable models

The hypothesis of error invariance is central to the instrumental variable literature. It means that the error term of the model is the same across all potential outcomes. In other words, this assumption signifies that treatment effects are constant across all subjects. It allows to interpret instrumental variable estimates as average treatment effects over the whole population of the study. When this assumption does not hold, the bias of instrumental variable estimators can be larger than that of naive estimators ignoring endogeneity. This paper develops two tests for the assumption of error invariance when the treatment is endogenous, an instrumental variable is available and the model is separable. The first test assumes that the potential outcomes are linear in the regressors and is computationally simple. The second test is nonparametric and relies on Tikhonov regularization. The treatment can be either discrete or continuous. We show that the tests have asymptotically correct level and asymptotic power equal to one against a range of alternatives. Simulations demonstrate that the proposed tests attain excellent finite sample performances. The methodology is also applied to the evaluation of returns to schooling and the effect of price on demand in a fish market.

Introduction
We are interested in the effect of a (possibly continuous) random treatment Z (with support Z ⊂ R p ) on a scalar outcome outcome Y . Let Y (z) be the scalar potential outcome of the outcome Y under treatment status z ∈ Z. The potential outcomes are not observed, but Y = Y (Z) is part of the data. We write where U (z) is an error term such that E[U (z)] = 0. Letting U = U (Z), we have Y = ϕ(Z) + U . This is a separable model (as opposed to nonseperable quantile models, see e.g. Chernozhukov and Hansen (2005)). The goal is to identify the causal effect ϕ, which allows to obtain the average treatment effect of changing treatment level from z ∈ Z to z ∈ Z, We say that Z is endogenous when the causal effect ϕ is not characterized by the distribution of Y conditional on Z. Even in this setting, it is possible to identifiy ϕ thanks to an instrumental variable W with support W ⊂ R q . Loosely speaking, an instrumental variable satisfies three conditions. First, it is (in some sense specified in the paper) independent of {U (z)} z . Second, it affects (causally) Y only through Z. Third, it is sufficiently related to Z.
Our goal is to test the following hypothesis H 0 : U (z) = U for all z ∈ Z.
We call this assumption "error invariance" since it stipulates that the error term does not vary with the treatment level.
There are two reasons why one would be interested in testing H 0 . First, under the instrumental variable conditions, several works in the instrumental variable literature (Newey and Powell (2003a); Horowitz and Lee (2007); Darolles et al. (2011); Chen and Pouzo (2012) among others) show that the causal effect ϕ is identified if H 0 holds. When H 0 is not satisfied, instrumental variable estimators can be more biased than naive estimators ignoring confounding since the error term U = U (Z) of the model can be dependent of W (through Z). Second, if H 0 is true, the treatment effects are constant (that is Y (z ) − Y (z) is a deterministic variable). This allows to estimate all counterfactuals for each observation in the dataset (individualized predictions) using an estimator of ϕ. It also means that the treatment does not distort the ranks of potential outcomes, so that the inequalities between subjects are not modified by the treatment.
The contributions of the paper are as follows. We propose two tests for H 0 . The tests rely on a scalar covariate X, which satisfies certain conditions outlined in the paper. The first test relies on the assumption that the potential outcomes are linear in the regressors.
In a first step, it estimates U as the residual of a two stage least-squares regression. Then, it tests that the estimator of U is uncorrelated with X(W k − E[W k ]) (for some component W k of W ), which is an implication of rank invariance under our assumptions. The second test is nonparametric. In this case, the error term U is estimated by the residual of a Tikhonov regression. The second step of the test is then the same as that of the linear test. In both cases, the treatment can be discrete or continuous. We study the asymptotic distribution of the test statistic and the power of the test. We discuss how to choose X: to maximize the likelihood that the conditions on X are satisfied, we argue in the paper that X should be chosen so as to be independent of W . Finally, the empirical performance of the tests is assessed through simulations and illustrated thanks to applications on returns to schooling and price elasticity of demand in a fish market.
Related literature. Our paper is related to a large literature on testing for (various implications of) constant treatment effects: see Koenker and Xiao (2002); Crump et al. (2008); Chernozhukov and Fernández-Val (2005); Ding et al. (2016); Hsu (2017) ;Goldman and Kaplan (2018); Chung and Olivares (2021); Sant'Anna (2021); Dai and Stern (2022) among others. Apart from Sant'Anna (2021), none of these papers considers the case where there is endogeneity. Sant'Anna (2021) only studies the setting where both the treatment and the instrument are binary and a monotonicity assumption such as in the literature on local average treatment effects holds (see Angrist et al. (1996)). Instead, our paper allows for confounding and continuous treatments.
Another related literature is that on the hypothesis of rank invariance (see Chernozhukov and Hansen (2005)). This assumption is the counterpart of our error invariance condition in nonseparable instrumental variable models (in contrast with the present paper which studies separable models). Chernozhukov and Hansen (2005) notes that identifica-tion of ϕ can be obtained under a weaker (but less interpretable) form of rank invariance called rank similarity. We could define an assumption analogous to rank similarity in our case of separable models. Our tests would also work for this assumption. We chose to focus solely on error invariance to simplify the exposition. Moreover, Dong and Shen (2018) and Frandsen and Lefgren (2018) develop tests of the rank invariance assumption in nonseparable instrumental variable models in the case where the treatment and the instrument are binary. The case with continuous treatment is much more challenging because regularization is needed for nonparametric estimation. Note also that Dong and Shen (2018) relies on the usual monotonicity assumption from local average treatment effects literature (see Angrist et al. (1996)), while our approach does not require it. The test developed in Frandsen and Lefgren (2018) bears similarities with our test since they are based on equivalent single restrictions. However, our model and the first step of our test are different. Moreover, Frandsen and Lefgren (2018) lacks a coherent theoretical framework justifying the test.
For instance, no assumption on the covariable (X in the present paper) is formulated in Frandsen and Lefgren (2018), so that it is unclear under which conditions rank invariance implies the linear hypothesis tested by Frandsen and Lefgren (2018) and when the test has power. In fact, one can show that the approach of Frandsen and Lefgren (2018) is implicitly based on moment restrictions constraining the data generating process in a similar way as our Assumption 2.1 on X, but Frandsen and Lefgren (2018) does not mention these hypotheses. In contrast, the present paper proposes a more formal framework in which it is possible to understand when the test will or will not work.
Outline. Section 2 presents the linear test along with theory, simulations and an application to returns to schooling. Then, in Section 3, we introduce the nonparametric test, discuss its asymptotic properties, evaluate its finite sample performance through numerical experiments and illustrate the test with an application to demand estimation. All the proofs along with some examples regarding the analysis of the power of the nonparametric test can be found in the supplementary material.
2 Linear test

Model
In this section, we develop theory for a test assuming that the mapping ϕ is linear and the for all z ∈ Z. We also impose the usual relevance assumption for instrumental variables in the linear model, that is E[W Z ] has rank equal to p. The aim is to recover β, since the average treatment effect of changing treatment level from z ∈ Z to z ∈ Z can be expressed To test for H 0 , we leverage a scalar covariate X satisfying the following condition: Assumption 2.1 There exists k ∈ {1, . . . , q} such that, for all z ∈ Z, This condition means that one component W k of the instrument is uncorrelated with the product of the unobserved heterogeneity of the model U (z) and the covariable X.
The tests of H 0 that we propose depend crucially on Assumption 2.1. Contrarily to H 0 , hypothesis 2.1 can be heuristically justified. Indeed, Assumption 2.1 is implied by the stronger condition where the sign ⊥ ⊥ stands for statistical independence. Hence, the econometrician should try to select a covariable X which is likely to be (jointly with {U (z)} z ) independent of W k .
The fact that {U (z)} z is independent of W corresponds to what is expected from a "good" instrument. As a result, in applications, we recommend to pick a variable X independent of W k since this should make Assumption 2.1 more likely to hold (the joint independence in (1) is not much stronger than the marginal independence, {U (z)} z ⊥ ⊥ W and X ⊥ ⊥ W k ). Finally, remark also that X does not need to be an instrument. In particular, its dependence with {U (z)} z and Z is left unrestricted.
Let us now define the population analog of the two-stage least squares estimator (henceforth, TSLS) in this context: In this case, under H 0 , β can be estimated by TSLS, that is β T SLS = β. However, when H 0 does not hold, then the bias of the TSLS estimator can be larger than that of the ordinary least squares (henceforth, OLS) estimator as illustrated by the following example.

Testable implication
This is the testable implication of H 0 that we use to build our test.
Remark that the role of the covariable X appears clearly in (2). Indeed, when W contains an intercept, the equality always holds by definition of U T SLS , regardless of the validity of H 0 . Hence, a test based on (3) would have no power. The role of X is therefore to give power to our test.

Sample test
Let us formally outline the test. Consider an i.i.d. sample {Y i , Z i , W i , X i } n i=1 generated from the model of Section 2.1. Let β T SLS be the TSLS estimator in this sample, given by The test statistic is then We use bootstrap to obtain the p-value of the test, since the variance of the statistic has a tedious expression and variance estimators based on analytical formulas tend to perform poorly in the presence of heteroscedasticity (see e.g. MacKinnon and White (1985)). The procedure to compute the p-value is as follows.
1. Draw n observations with replacement from the sample 2. Compute the bootstrapped statistic T * n,b on the bootstrapped sample 3. Repeat steps 1-2 B times (with B large) so as to get the collection of bootstrapped If the p-value so obtained is smaller than α (the nominal size of the test), the null hypothesis is rejected at the α nominal level. Notice that in the above procedure we are computing a symmetric p-value. An alternative procedure is to compute an equal-tailed p-value, but in our simulations the symmetric one has a satisfying performance.

Power analysis
The test will have asymptotic power equal to 1 under alternative hypotheses for which Since We see that there are two reasons why the test would have power: (1) the variable W k is correlated with U (Z)X; (2) it holds that β = β T SLS and E[X(W k −E[W k ])Z ] = 0 for at least one in {1, . . . , p}.
We argue that (1) and (2)  Overall, although there are cases where the test does not have power, the above discussion suggests that the test has power against a wide range of alternative hypotheses. The following example illustrates this claim.
Example 2. Let (W, E, X) be a 3 × 1 Gaussian vector with mean zero and variance equal to the identity matrix. Let also Z = W + W E + ρW X, U (Z) = Z(E + ρX), where ρ ∈ R, and β = 0 so that Y = U (Z) = Z(E + ρX). In this case, we have In the present context, it is well- As a result, the test has no power only when X and U are independent, that is ρ = 0, which is a degenerate case.
In this parametric example, we see that the alternatives under which the test has no power have measure equal to 0 for the uniform measure. The example also shows the role of X in providing power to the test, since the latter does not have power only when X and (Z, {U (z)} z ) are independent.

Asymptotic theory
In this section, we state the asymptotic properties of the test statistic. We make the following assumption, which ensures the convergence of the TSLS estimator.
exists and has rank p and E[W W ] exists and has full rank.
We have the following theorem.
Theorem 2.1 Let Assumptions 2.1 and 2.2 hold. Then, T n converges to a zero mean Gaussian distribution under H 0 , while |T n |/ √ n → C = 0 when (4) holds.
The influence function representation of the test statistic is given in Lemma 1 in the supplementary material. The asymptotic variance of the test statistic can be derived from this result. Note that the fact that |T n |/ √ n → C = 0 when (4) holds implies that the power of the test goes to 1 as n goes to ∞ under alternative hypotheses satisfying (4).

Simulations
We study the following data generating process. Let (W 2 , E, X, Z 3 ) follow a standard 4-dimensional Gaussian distribution. We define The variable Z 2 suffers from confounding and is instrumented by where ρ ∈ R. The null hypothesis H 0 holds when ρ = 0. The outcome is Y = U (Z) so that the causal regression function of interest is ϕ = 0. We set k = 2, that is we use the second component of W to compute the test statistic We generate data with sample sizes n ∈ {100, 200, 500, 1000}. We investigate the empirical size of the test when ρ = 0 ( Table 1). The results are averages over 10000 replications using B = 1000 bootstrap resamples. We also study the empirical power of the test when ρ = −1, −0.9, . . . , 1, with a thousand replications and bootstrap resamples (See Figure 1).
The empirical size of the test is almost nominal even for low sample sizes. The power of the test increases as the deviation from the null (|ρ|), or the sample size, become larger. n 100 250 500 1000 Empirical size at 5% 0.0397 0.0512 0.0506 0.0506 Empirical size at 10% 0.1038 0.1050 0.1038 0.1014 Table 1: Empirical size of the test of theoretical size 5% and 10% for various sample sizes.

Application on returns to schooling
In this section, we test for error invariance in a famous application on returns to schooling from Card (1995). The database is available in the R package ivreg and is an extract from the 1976 National Longitudinal Survey (NLS) of young men. In this example Y = log(wage), Z = (1, education, experience, experience 2 , smsa, south, ethnicity) , where 1 denotes a constant, education and work experience are measured in years and smsa, south, ethnicity are controls whose definition can be found in the ivreg package. Education and experience are endogenous since they depend on individual's ability which is unobserved. The instrument is W = (1, nearcollege, age, age 2 , smsa, south, ethnicity) , where nearcollege is an indicator whose value is equal to 1 when the individual grew up near an accredited four-year college. The validity of the instrument is documented in Card (1995). This corresponds to one of the specifications in the original paper of Card (1995). The variable X is an indicator for being married. It is unlikely that being married is (jointly with the unobserved heterogeneity of the model) dependent from growing up near a four-year college. In fact, a χ 2 test does not reject hypothesis of independence of W 2 = nearcollege and X = married (p-value equal to 0.55). Hence, we can reasonably

Nonparametric test
In this section we extend the linear framework to consider a nonparametric treatment function. We keep the notation previously introduced and write The error term is now supposed to be mean independent of the instrument ( Let L 2 (Z) be the set of functions that are square integrable with respect to the distribution of Z. We maintain the following assumption which is also standard in the NPIV literature.
Assumption 3.1 The following holds: The first part of the assumption is the completeness condition standard in the NPIV literature. It means that W is sufficiently associated with Z. Under the null hypothesis H 0 , such a condition is necessary and sufficient for identifying ϕ (see, among others, Carrasco et al. (2007), Darolles et al. (2011), Newey andPowell (2003b), D'Haultfoeuille (2011)).
The second part states that the NPIV model in the above equation is well specified. We will discuss this assumption in Section 3.2.
Recall that the hypothesis that we want to test is Under the null hypothesis of error invariance the treatment function ϕ must satisfy the following equation Under the completeness condition (Assumption 3.1 (i)) the solution to the above equation is unique and can be consistently estimated (we introduce an estimator below). Thus, under the null hypothesis of error invariance we can consistently estimate the nonparametric treatment function ϕ. Differently, when H 0 does not hold, the solution to Equation (6) will be different, in general, from the treatment function of interest ϕ. This solution, denoted ϕ, is introduced in Assumption 3.1 (ii). So, if H 0 is not true, although we will be able to estimate consistently the solution to Equation (6), we will in general not obtain a consistent estimator of ϕ. It is therefore crucial to test for the error invariance hypothesis.
Similarly to the previous section, the testable implication of the error invariance assumption is and ϕ denotes the unique solution to Equation (6)  Let ϕ be an estimator for ϕ (below we define our estimator), and let us define the esti- Our statistic is the empirical counterpart of the above moment In the following section, we explain how to compute the above statistic and implement the test based on it. In Section 3.2, we present the assumptions and the results for the validity of the test. Next, in Section 3.3, we provide evidence about the finite sample performance of the test based on S n . Finally, we illustrate the procedure by an application to demand estimation in Section 3.4.

Construction of the statistic and implementation of the test
In this section, we outline the computation of the statistic S n and implement the test. To this end, we need (i) to compute ϕ(Z i ) for i = 1, . . . , n, and (ii) to obtain the p-value necessary for testing.
Computation of ϕ. Let π and τ be positive functions on R. We denote with L 2 π (R) the space of square integrable functions with respect to π, and let L 2 τ (R) be similarly defined. The functions π and τ are introduced for the aim of generality and for technical reasons.
In simulations they are set to Gaussian density functions. Let us briefly explain their roles.
First, from a technical point of view it would be ideal to work with the spaces L 2 (Z) and L 2 (W ), i.e. the spaces of square integrable functions with respect to the distribution of Z and W . However, since we do not know such distributions, it is convenient to replace L 2 (Z) and L 2 (W ) with L 2 τ (R) and L 2 π (R) and work with the latter spaces. Second, by working with L 2 τ (R) and L 2 π (R) we can obtain that the estimator of the dual of A (defined below) is actually the dual of the estimator of A. This is an important property for establishing the asymptotic normality of our statistic S n under H 0 . Let f W stand for the density of W with respect to the Lebesgue measure, and assume that τ (w) > 0 whenever f W (w) > 0.
Then, Equation (6) To build an estimator for ϕ we start from the above integral equation. Let us denote by f ZW the joint density of (Z, W ) with respect to the Lebesgue measure and introduce L 2 π⊗τ (R 2 ), which is the set of functions from R 2 to R square integrable with respect to the product measure π ⊗ τ (the measure generated by the product density π · τ ). By assuming .

Its Hilbert adjoint
We also define the "left hand side" of the integral equation (7) as We estimate A and A * by replacing f ZW with its kernel estimator where K Z and K W are two kernel functions while (h Z , h W ) denote two bandwidths converging to zero as the sample size increases. The mapping r is instead estimated by To estimate ϕ, i.e. the solution to Equation (7), we use a Tikhonov regularization scheme where I stands for the identity operator, while λ > 0 denotes the Tikhonov regularization parameter. See, e.g.,  and Darolles et al. (2011).
To compute ϕ we need to select two bandwidths (h Z , h W ) and a regularization parameter λ. We will explain below how to do it. For the moment, let us outline the computation of ϕ in the previous display for given (h Z , h W , λ). Although ϕ seems to have an abstract expression, its computation is straightforward as it reduces to matrix products. To see this, we first approximate A and A * by using a common bias computation similarly as in Centorrino et al. (2017) : The above approximations hold for h Z , h W → 0. Let M Z be the n × n matrix having on the ith row and jth column the element K((Z i − Z j )/h Z )/[π(Z i )nh Z ] and let M W be the n × n matrix having on the ith row and jth column the element K((W i − . Let − → r and − → Aφ (for a generic function φ) be similarly defined. By the expression of ϕ and the above approximations we get up to an approximation error. Thus, the computation of ϕ at the data points reduces to a simple matrix computation and the residuals can be easily computed as Let us turn now to the selection of the smoothing parameters (h Z , h W ) and the regularization parameter λ. The bandwidths h Z and h W are selected by the Silverman Rule of Thumb, so that h Z = σ Z n −1/5 and h W = σ W n −1/5 , where σ Z and σ W denote the sample standard deviations of Z and W . Finally, to select the Tikhonov regularization parameter, we use the Cross-Validation method developed in Centorrino et al. (2017). Thus, where Here A −i and A * −i denote the "leaveone-out" versions of A and A * that use the entire sample except the ith observation. By the approximations of A and A * previously outlined, the vector where diag(M W ) denotes the diagonal matrix having the same main diagonal as M W , and diag(M Z ) is similarly defined. The objective function in (9) has a U-shaped form as a function of λ and can be minimized in a simple way by numerical methods, see Centorrino et al. (2017).
For the sake of clarity, we sum up the steps needed for computing S n as follows: 1. Select λ according to the Cross-validation method in Equation (9) 2. Use λ to compute { ϕ(Z i ) : i = 1, . . . , n} as in (8) and compute the residuals as 3. Compute the statistic as Implementation of the test. To implement the test, we are just left with the computation of the p-value. As we show in the next section, under the null hypothesis the statistic S n is asymptotically normal but the asymptotic covariance has an intricate expression.
So, to obtain the p-value necessary for testing we rely on the pairwise bootstrap. The validity of the pairwise bootstrap is confirmed by our simulations and is not surprising given that we show asymptotic normality. In fact, under different conditions, Chen et al. (2003) show the validity of the bootstrap for a a real-valued estimator based on a first step infinite-dimensional estimator, as is ϕ in this paper.
The procedure goes as follows: 1. Draw n observations with replacement from the sample 2. Compute the bootstrapped statistic S * n,b on the bootstrapped sample using the same bandwidths (h Z , h W ) and regularization parameter λ as in the original sample 3. Repeat steps 1-2 B times (with B large) to obtain the collection of bootstrapped

Asymptotic behavior
In this section, we outline the assumptions under which we obtain the asymptotic properties of the test statistics. The following definition introduces some regularity features that the joint density f ZW and the NPIV function ϕ must satisfy.
Definition 1 For a given function γ and for α ≥ 0 , s > 0, the space B s,α γ (R ) is the class of functions g : R → R satisfying: g is everywhere (m − 1) -times partially differentiable for m − 1 < s ≤ m; for some R > 0 and for all x, the inequality holds true, where Q = 0 when m = 1 and Q is an (m − 1)th degree homogeneous polynomial in (y − x) with coefficients the partial derivatives of g at x of orders 1 through m − 1 when m > 1; ψ is uniformly bounded by a constant when α = 0 and the functions g and ψ have finite αth moments with respect to 1/γ when α > 0, i.e. |g(x)| α /γ(x) dx < ∞ and Let f Z denote the density of Z with respect to the Lebesgue measure and let us define Assumption 3.2 ϕ ∈ B ρ,0 1 (R) ∩ L 2 π (R) and f ZW /(π τ ) ∈ B ρ,0 Assumption 3.2 is a common regularity condition allowing for several degrees of integrability and differentiability . Under the above assumption A : L 2 π (R) → L 2 τ (R), A * : L 2 τ (R) → L 2 π (R), and r ∈ L 2 τ (R). Notice that A * is actually the Hilbert adjoint of A. This aspect is used multiple times in the proofs. Assumption 3.3 K Z and K W are symmetric kernels of order ρ with bounded support.
The kernel orders in Assumption 3.3 are assumed to be equal only for notational simplicity.
To simplify our theoretical exposition we further assume that h W = h Z and denote each bandwidth by h. However, our proofs also hold when such smoothing parameters are set to different values (and the kernel have different orders).
The order ρ of the kernels K Z and K W , the bandwidth h and the regularization parameter λ satisfy the following assumption.
Such conditions allow us to control the regularization bias and the variance of ϕ appearing in the expansion of S n . In particular, nh 2 λ 3/2 allows controlling the "variance term". nλ 2 and h ρ λ −3/4 allow us to obtain the negligibility of the regularization bias that appears in the expansion of S n and is due to the ill-posedness of the inverse problem.
Let us denote with ·, · the inner product of either L 2 π (R) or L 2 τ (R), and let · be the norm induced by such an inner product. The specific inner product or norm we refer to will be clear at each time from the context. Let (λ j , ϕ j , ψ j ) j be the Singular Value decomposition of the operator A, where (λ j ) j is the sequence of singular values in R, (ϕ j ) j is an orthonormal sequence in L 2 π (R), and (ψ j ) j is an orthonormal sequence in L 2 τ (R). The following assumptions introduces the usual source conditions. Assumption 3.5 Let ·, · denote the inner product on L 2 π (R).
Source conditions are standard in the NPIV literature (see  or Darolles et al. (2011)). In general, the source condition is imposed only on ϕ when the interest is on estimating ϕ. Here, we need a source condition on both ϕ and g to establish the asymptotic properties of the statistic based on the residuals from the nonparametric regression. To this end, we also need some regularity conditions on the estimator ϕ. So, let N [·] ( , Φ, || · || P ) be the bracketing number of size of a class of functions Φ, where || · || P denotes the L 2 norm with respect to the probability P, see .
Assumption 3.6 There exists a class of functions Φ such that The above condition is a high-level assumptions allowing us to handle an empirical process in the expansion of our statistic. It essentially introduces a Donsker feature for the function ϕ. Such an assumption has been used among others by Rothe (2009)   or Van der Vaart and Wellner (1996). We have also derived alternative proofs based on sample splitting or cross-fitting that avoid such a high-level condition. We have run simulations with cross-fitting and sample-splitting, but we did not notice any improvement in terms of finite sample performance. Thus, we choose to keep the above assumption and avoid sample splitting (or cross-fitting) for a simple implementation of the test. Alternatively such a high-level condition could be avoided by using a Sobolev penalized method for estimating ϕ. This however would produce much longer proofs and complicate the implementation of the test.
Let P n be the empirical mean operator. We have the following theorem.
Theorem 3.1 Let Assumptions 2.1 and 3.1 to 3.6 hold. Under H 0 , The proof of the Theorem is in Section B.1 of the supplementary material. Let us briefly comment on the asymptotic expansion (influence-function representation) of the statistic.
The first term on the right hand side is the version of our statistic that we could use had we observed the error U and E[W ]. The second term arises because the expectation E[W ] is unobserved and is replaced by its estimated counterpart. Similarly, the third term arises because the error U is unobserved and to estimate it we replace the true function ϕ 0 with its estimator ϕ. More in detail, such a term originates from a (nonparametric) bias involving the difference ϕ − ϕ that appears in the expansion of S n . This inflates the asymptotic variance of the statistic, as compared to the case where U is known, and hence represents the price to pay for not observing the error term. To handle this bias term, we rely on decompositions from Darolles et al. (2011) and Babii and Florens (2017). The expansion in Theorem 3.1 allows us to obtain the asymptotic normality of our statistic under the null hypothesis. To this end, notice that the first term has zero expectation under the null of error invariance. The second term has clearly a null expectation. Finally, since E[U |W ] = 0, the last term also has zero expectation. Therefore, a standard Central Limit Theorem implies the asymptotic normality of S n . As for the linear test, the fact that
To check the robustness of our test with respect to the choices of the smoothing and regularization parameters, we let them vary around the benchmark values. So, we set are the bandwidths set according to the Silverman's Rule of Thumb while λ * is the regularization parameter selected according to the Cross-Validation method, as seen in Section 3.1. C h and C λ are fixed constants.
We run simulations for C h , C λ ∈ {0.5, 1, 2} and we consider the modest sample sizes of n = 100, 250, 500. As a kernel we use the standard Gaussian density. We set π to the Gaussian density with mean equal to the sample mean of Z and variance equal to twice the sample variance of Z. Similarly, the measure τ is set to the Gaussian density with mean equal to the sample mean of W and variance equal to twice the sample variance of W .  used a similar strategy to select π and τ in their simulation setting. Intuitively, since π and τ have to be two densities and at the same time appear at the denominators in A and A * , we want them (i) to integrate to 1 and (ii) to converge towards zero sufficiently slowly on the tails. To speed up computations, we use the warpspeed method by Giacomini et al. (2013): for each Monte Carlo iteration we draw a single bootstrap sample, and we use the whole set of bootstrap statistics to compute the bootstrap p-values associated with each original statistic. We perform a very large number of Monte Carlos iterations equal to 10,000.
The results under the null hypothesis of error invariance (γ = 0) are reported in Table   2. The tests are implemented at 5 and 10 percent nominal levels. The result show that the test is reasonably stable with respect to the choices of λ, h Z , and h W and controls well the size under the null hypothesis. Also, as expected, the error in the rejection probability To sum up, these simulations show that (i) for the moderate sample sizes of n = 100, 250, 500 the tests display a good performance both in terms of size control and in terms of power, and (ii) the results are stable with respect to the choices of the smoothing and regularization parameters.

Application to real data
In this section we apply the test of error invariance to a fish demand equation. This describes the relationship between the demanded quantity of fish and its price. It is common practice in econometrics to assume that prices are endogenous and estimate demand equations by instrumental variable regressions.
The data we use come from Graddy (2006) and contain daily observations at the New York Fulton fish market about the log of the price, the sold quantity, an indicator for each day of the week, and the wind speed registered in the previous day. The dataset can be downloaded at "https://users.ox.ac.uk/ nuff0078/EconometricModeling/". The total number of observations is 111.
We study a demand equation. The outcome Y is equal to the logarithm of the quantity sold and Z is the logarithm of the market price. Market price is likely to be endogenous since it also depends on expected demand, which itself also affects the quantity. Following the literature (Graddy (2006)), we choose a variable related to weather as an instrument, so we set W equal to the wind speed recorded on the day corresponding to the observation.
Such a variable is viewed as sufficiently correlated with the price (Z), since the weather affects the ability to fish, and is at the same time exogenous with respect to the errors {U (z)} z since it is probably unrelated to demand shocks related to human factors.
The choice of the covariable X is crucial. We set X equal to the indicator for week days (Monday to Thursday). Such a variable should not be correlated with the wind speed since the weather does not depend on the day of the week. In fact, a Kolmogorov-Smirnov test does not reject the hypothesis of independence between X and W (p-value of 0.17).
Hence, we can reasonably assume that E[U (z)X(W − E[W ])] = 0 for all z ∈ Z, which corresponds to Assumption 2.1. Table 3 reports the p-values of the test. Similarly as in Section 3.3, to check the robustness of our findings we perturb the benchmark bandwidths (h * W , h * Z ) and the benchmark regularization parameter λ * . We use 1, 000 bootstrap resamplings. Our test rejects the null hypothesis of error invariance at the 5 percent nominal level for most configurations considered. So, we conclude that the causal intepretation of the separable NPIV regression should be doubted for this specific example.
It holds that Hence, (A.1) is a linear projection of Y on Z. Since this linear projection is unique (because

A.2 Proof of Theorem 1
Theorem 1 is a direct corollary of the following lemma, which gives the influence function representation of the test statistic.
Lemma A.1 Let Assumptions 2.1 and 2.2 in the main text hold. Then, we have We deal with the first term on the right-hand side of (A.2). We have Then, we handle the second term on the right-hand side of (A.2). It holds that Since β T SLS − β = O P (1/ √ n), and by the law of large numbers Moreover, by standard computations, we have which leads to the result.

2
B On the nonparametric test B.1 Proof of the results of Section 3 Recall that P n represents the empirical mean operator. Let us denote with P the population mean operator, so Pf (ξ) = Ef (ξ) for any random variable ξ. We finally denote with · π the norm of the L 2 π (R) space, i.e. f 2 π = |f (z)| 2 π(d z) for any f ∈ L 2 π (R).

B.1.1 Proof of Theorem 2
Let us decompose the statistic as follows: We start by showing that the last term on the RHS equals E[U X] · √ nP n (W − E[W ]) + o P (1). Notice that √ nP n U X(W − E[W ]) = P n U X · √ n(P n − P)W . By Assumption 3.6 in the main text, with probability approaching one It follows from the definition of bracketing . Assumption 3.6 ensures that the RHS of this inequality is finite for any . So, the class Υ is Glivenko-Cantelli and by Theorem 19.4 in  we get that (P n − P)φ = o P (1) uniformly in φ ∈ Υ. This and the previous display lead to P n U X = P U X + P (ϕ − ϕ) X + o P (1) .
To show the negligibility of the second term on the RHS, by the Law of Iterated Expectations and the Cauchy-Schwartz inequality we get where · π denotes the norm on L 2 π (R). The second factor on the RHS is finite by Assumption 3.2 in the main text. The first factor is instead o P (1) by Lemma B.5. By putting together these results, we find that We now show the negligibility of the second term on the RHS of Equation (B.1). Combining the Law of Iterated Expectations and Assumption 3.2 leads to with ϕ−ϕ 2 π = o P (1), as previously found. So, by Assumption 3.6 and the boundedness of To show the behavior of the statistic under the alternative, we can proceed as at the beginning of this proof to obtain where U = Y − ϕ(Z) and the first leading term is different from 0 under H 1 .

B.1.2 Auxiliary lemma
In this section we provide an auxiliary lemma we used in the proofs of Theorem 2. Let us denote with ·, · the inner product on either L 2 π (R) or L 2 τ (R), and let · denote the norm induced by such an inner product. Whether ·, · and · refer to either L 2 π (R) or L 2 τ (R) will be clear by their arguments. Also, given the operators A and A * we define Aφ 2 and A * 2 op := sup Since A and A * are linear bounded operators, both A 2 op and A * 2 op are finite. Also, since A * is the Hilbert adjoint of A, it holds that A 2 op and A * 2 op are equal.
Lemma B.1 Under Assumptions 3.1 to 3.6 in the main text, we have Proof. Let us start with the following decomposition The previous decomposition holds because where the second line follows from the Cauchy-Schwartz inequality and the continuity of A−A, while the third line follows from Assumption 3.5(i) and Lemma B.4(iv). Assumption 3.5(i) ensures that η ≥ 2. As already noticed before, || A − A|| op / √ λ = o P (1) and || r − A ϕ||/ √ λ = o P (1). So, by the above display and Assumption 3.4 we find We now handle the term √ n ϕ λ − ϕ, g . By Assumption 3.5(ii) and Lemma B.
By proceeding along the same lines we find that E 4th term RHS of (B.17) 2 ≤ Ch 4ρ n = o(1) .
We are now able to obtain the Influence function representation for √ n Ξ 1 , g in Equation (B.16). By the previous five displays and since (A * A) −1 g < ∞, the first term on the RHS of (B.16) equals To conclude the proof, notice first that Equation (B.17) and the four displays after it ensure that √ nA * ( r − A ϕ) = O P (1). So, the norm of the second term on the RHS of (B.16) is upperbounded as follows By recalling that we have set η = 2, we have (A * A) η/2 = A * A. Thus, by definition of f we get g = (A * A) η/2 f = A * A f and Since (λI + A * A) −1 A * is a regularization scheme, the RHS of the above display converges to 0 as λ → 0 (see Kress (1999)). This concludes the proof.

B.1.3 Auxiliary Results
The following lemma is borrowed from  (see their Lemma A.1(b)).
Lemma B.2 Under Assumptions 3.2 to 3.4 in the main text, we have The following lemma is well established in the literature on empirical process theory. A proof can be found in Andrews (1994). The lemma that follows gathers several useful results about inequalities of norms involving compact operators. Its proof can be found in Florens et al. (2011) and.