module_for_agricultural_supply_at_regional_level
Differences
This shows you the differences between two versions of the page.
Both sides previous revisionPrevious revisionNext revision | Previous revision | ||
module_for_agricultural_supply_at_regional_level [2020/03/11 10:46] – [Detailed discussion of the equations in the supply model] matsz | module_for_agricultural_supply_at_regional_level [2023/09/08 12:11] (current) – [Price depending crop yields and input coefficients] massfeller | ||
---|---|---|---|
Line 89: | Line 89: | ||
Where “asym” is the land asymptote, i.e. the maximal amount of economically usable agricultural area in a region when the agricultural land rent goes towards infinity. For an application where the land market is used see Renwick et al. (2013). | Where “asym” is the land asymptote, i.e. the maximal amount of economically usable agricultural area in a region when the agricultural land rent goes towards infinity. For an application where the land market is used see Renwick et al. (2013). | ||
- | Set aside policies have changed frequently during CAP reforms. The recent specification is covered in the context of the premium modelling in Section | + | Set aside policies have changed frequently during CAP reforms. The recent specification is covered in the context of the premium modelling in Section |
\begin{align} | \begin{align} | ||
Line 130: | Line 130: | ||
& \forall r,n,j | & \forall r,n,j | ||
\end{split} | \end{split} | ||
- | \end{align} | + | \end{align} |
Indices: | Indices: | ||
+ | |||
+ | \(r\) = region \\ | ||
+ | \(i\) = crop \\ | ||
+ | \(j\) = crop group \\ | ||
+ | \(k\) = technological crop option (high/low yield) \\ | ||
+ | \(n\) = nutrient (N/P/K) \\ | ||
+ | \(isPerm_j\) = indicates that crop group \(j\) contains permanent crops | ||
+ | |||
+ | Endogenous choice variables: | ||
+ | |||
+ | \(levl_{rik}\) = Area (ha) of each crop \(i\) and technology \(k\) in region \(r\). | ||
+ | \(fmine_{rnj}\) = Application of mineral fertilizer \(n\) to crop group \(j\) in region \(r\). | ||
+ | \(fexcr_{rnj}\) = Application of manure \(n\) to crop group \(j\) in region \(r\). | ||
+ | \(fcrex_{rnj}\) = Allocation of crop residue \(n\) to crop group \(j\) in region \(r\). | ||
+ | |||
+ | Parameters: | ||
+ | \(ret_{rni}\) = Retention (uptake) of nutrients by the crop \\ | ||
+ | \(res_{rni}\) = Crop residues output \\ | ||
+ | \(biofix_{rni}\) | ||
+ | \(λ_{rnik}^{prop} \) = Over-fertilization factor, calibrated \\ | ||
+ | \(λ_{rni}^{const} \) = Over-fertilization term, calibrated \\ | ||
+ | \(soil_{rn}\) = Soil factor \\ | ||
+ | \(yf_{rnik}\) = Yield factor for technologies \\ | ||
+ | \(loss_{rn}\) = Loss rate \\ | ||
+ | \(ϕ_{r", | ||
+ | \(ϕ_{r," | ||
+ | |||
+ | The reader may have noted that there is no loss rate for manure in the Equation 96 FIXME . CAPRI does contain such loss rates, but they are specific for each animal type and therefore happens on the manure supply side of the regional manure balance (see section on input allocation). | ||
+ | |||
+ | The model contains three types of manure: N-manure, P-manure and K-manure. From an agricultural point of view this may seem odd. It might be more intuitive to think of one type of manure per animal category. The motivation is to keep the system simple and flexible. With the present representation, | ||
+ | |||
+ | The supply of each manure type is collected in a “pool” for each regional farm model, i.e. for each NUTS2 region. Regions within a member state may trade manure, subject to a cost. The supply in the pool plus the traded quantities has to be distributed to the crops in the region, i.e. there is an equality-restriction in place. This is handled in the equations “FertDistExcr_” and “ManureNPK_”. Note that fertilizer flows are measured in //tons//, for the sake of scaling, whereas other total quantities in CAPRI are measured in 1000 tons. Hence the factors 1000 and 0.001. | ||
+ | |||
+ | \begin{equation} | ||
+ | \sum_j fexcr_{rnj}= 1000 v\_ManureNPK_{rn} | ||
+ | \end{equation} | ||
+ | |||
+ | \begin{equation} | ||
+ | v\_ManureNPK_{rn} + \sum_s T_{rs}nutshr_{rn} = 0.001 \sum_{i\in Anim_j,k} levl_{rik}o_{rnik}(1-loss_{rin}) \quad \forall r, n | ||
+ | \end{equation} | ||
+ | |||
+ | where \\ | ||
+ | \(o_{rnik}\) is the output of manure nutrient \(n\) from animal type \(i\) using technology \(k\) in region \(r\), \\ | ||
+ | \(nutshr_{rn}\) is the average content of each nutrient in the regional manure pool, \\ | ||
+ | \(T_{rs}\) is the quantity of manure traded from \(r\) to \(s\), \\ | ||
+ | \(isAnim_{i}\) indicates that activity \(i\) is an animal production activity | ||
+ | |||
+ | Equation “FertDistMine_” allocates total mineral fertilizer sales to the crops / group of crops. | ||
+ | |||
+ | \begin{equation} | ||
+ | \sum_j fmine_{rnj} = -netPutQuant_{rn} | ||
+ | \end{equation} | ||
+ | |||
+ | Finally, crop residues and atmospheric definition are distributed in equation “FertDistCres_”. | ||
+ | |||
+ | \begin{equation} | ||
+ | \sum_j fcers_{rnj} = \sum_{i \notin isPerm_j,k} levl_{rik}res_{rni}(techf_{rink}+1) | ||
+ | \end{equation} | ||
+ | |||
+ | One flow from a source s={" | ||
+ | |||
+ | To resolve the ill-posedness of the fertilizer distribution, | ||
+ | |||
+ | To develop this probabilistic model, we assume that the decisions of the farmer are separable and taken in two steps: first, the farmer decides about the cropping plan and just ensures that the total amount of fertilizer available is sufficient. This is called the outer model. Then, a statistical model is solved that finds the most probable fertilizer flows out of the continuum of possible ones. This is called the inner model. The structure with outer and inner models makes the problem a bi-level programming one. | ||
+ | |||
+ | To implement the bi-level programming problem in a way that does not change the present structure of the model (with just one optimization solve of the representative farm model) we implement the inner model by its optimality conditions. By carefully choosing the proper probability density functions we ensure that no complementary slackness conditions are needed, so that the inner model is simple to solve. For this the gamma density function is very suitable, as it has a support from zero to infinity, with a probability that goes towards zero as the random variable goes to zero. | ||
+ | |||
+ | The parameters of the gamma function are determined in the calibration step, described further below, and then kept constant in simulation. The gamma density function for some random variable x has the form | ||
+ | |||
+ | \begin{equation} | ||
+ | p(x|\alpha, | ||
+ | \end{equation} | ||
+ | |||
+ | where Γ(α) is the gamma function, and α and β are parameters that determine the shape of the density function. The gamma density is nonlinear, and the joint density, being the product of the densities of all nutrient flows, is even more so. In order to reduce nonlinearity we note that are interested in finding the highest posterior density, i.e. maximizing a joint density function, and since the maximum is invariant to monotonous positive transformations we compute the logarithm of the joint density, which will be the sum of terms like the following (the constant term has been omitted since it also does not influence the optimal solution for x): | ||
+ | |||
+ | \begin{equation} | ||
+ | log \, p(x|\alpha, | ||
+ | \end{equation} | ||
+ | |||
+ | Maximization of the logged density under the constraints that the nutrient balance restrictions of the supply modes have to be met gives a set of equations that define explicit and unique fertilizer flows, where v are the Lagrange multipliers of the source-pool restrictions and u the Lagrange multipliers of the nutrient balance equation. | ||
+ | |||
+ | \begin{equation} | ||
+ | \text{FOC w.r.t. manure use:} \\ | ||
+ | \frac {\alpha_{r, | ||
+ | \end{equation} | ||
+ | |||
+ | \begin{equation} | ||
+ | \text{FOC w.r.t. crop residues use:} \\ | ||
+ | \frac {\alpha_{r, | ||
+ | \end{equation} | ||
+ | |||
+ | \begin{equation} | ||
+ | \text{FOC w.r.t. mineral feritilzer use:} \\ | ||
+ | \frac {\alpha_{r, | ||
+ | \end{equation} | ||
+ | |||
+ | The system of FOC contains expressions of the type \(1/fmine\) which is likely to impair performance as the second derivatives are not constant (CONOPT computes second derivatives). Therefore, the first term in each FOC was turned into a new variable \(z\) defined as \( z_{r," | ||
+ | |||
+ | ===Balancing equations for outputs=== | ||
+ | |||
+ | Outputs produced must be sold – if they are tradable across regions – or used internally, as in the case of young animals or feed. | ||
+ | |||
+ | \begin{equation} | ||
+ | \sum_{act}Levl_{r, | ||
+ | \end{equation} | ||
+ | |||
+ | In the case of quotas (milk, for sugar beet) the sales to the market may be bounded (noting that NETTRD = v_netPutQuant in the code): | ||
+ | |||
+ | {{:: | ||
+ | |||
+ | As described in the data base chapter, the concept of the EAA requires a distinction between young animals as inputs and outputs, where only the net trade is valued in the EAA on the output side. Consequently, | ||
+ | |||
+ | \begin{equation} | ||
+ | \sum_{aact}Levl_{r, | ||
+ | \end{equation} | ||
+ | |||
+ | In combination with the standard balancing equation shown above, the NETTRD variable for young animals on the output side becomes negative if the YANUSE variable for a certain type of young animals exceeds the production inside the region. | ||
+ | |||
+ | ===The objective function=== | ||
+ | |||
+ | The objective function is split up into the linear part, the one related to the quadratic cost function for activities, and the quadratic cost function related to the feed mix costs: | ||
+ | |||
+ | \begin{equation} | ||
+ | OBJE=\sum_r LINEAR_r+QUADRA_r+QUADRF_r | ||
+ | \end{equation} | ||
+ | |||
+ | The linear part comprises the revenues from sales and the costs of purchases, minus the costs of allocated inputs not explicitly covered by constraints (i.e. all inputs with the exemptions of fertilisers, | ||
+ | |||
+ | \begin{equation} | ||
+ | LINEAR_r= \sum_{io} NETTRD_{r, | ||
+ | \end{equation} | ||
+ | |||
+ | The quadratic cost function relating to feed is defined as follows: | ||
+ | |||
+ | \begin{equation} | ||
+ | QUDRAF_r = \sum_{aact, | ||
+ | \begin {matrix} | ||
+ | LEVL_{r, | ||
+ | (a_{r, | ||
+ | \end{matrix} | ||
+ | \right] | ||
+ | \end{equation} | ||
+ | |||
+ | The marginal feed costs per animal increase hence linearly with an increase in the feed input coefficients per animal. It should be emphasised that this is the main mechanism that “stabilises” the feed allocation by animals. The two balances on feed energy and protein alone would otherwise leave the feed allocation indeterminate and give a rather “jumpy” simulation behaviour. | ||
+ | |||
+ | There is another more complex PMP term (equation quadra_ in supply_model.gms, | ||
+ | |||
+ | A final term relates to the entitlements introduced with the 2003 Mid Term reviews. If those entitlements are overshot, a penalty term equal to the premium paid under the respective scheme (regional, historical etc.) is subtracted to the objective. Accordingly, | ||
+ | |||
+ | ===Sugar beet (M. Adenäuer, P. Witzke)=== | ||
+ | |||
+ | The Common Market Organisation (CMO) for sugar regulates European sugar beet supply with a system of production quotas, even after the significant reforms of 2006, up to year 2017 when the quota system expired. Before that reform, two different quotas had been established subject to different price guarantee (A and B quotas, qA and qB). Beet prices were depending on intervention prices and levies to finance the subsidised export of a part of the quota production to third countries. Sugar beets produced beyond those quotas (so called C beets) were sold as sugar on the world market at prevailing prices, i.e. formally without subsidies. However, a WTO panel initiated by Australia and Brazil concluded that the former sugar CMO involved a cross-subsidisation of C-sugar from quota sugar such that all exports of C sugar was also counted in terms of the EU’s limits on subsidised exports. As a consequence, | ||
+ | * A-beets receiving the highest price derived from high sugar prices (and before the 2006 reform less a small levy amount) | ||
+ | * B-beets receiving a lower price as the applicable levies were higher before the reform. However, the 2006 reform eliminated the distinction of A and B quotas. Furthermore, | ||
+ | * C-beets receiving the lowest price, formerly derived from world market sugar prices, now derived from ethanol prices. | ||
+ | |||
+ | The high price sector covers for farmers at least the farm level quota endowment. However, the sugar industry may grant high prices also for a limited, “desirable” over-quota production, for example to avoid bottlenecks in sugar or ethanol production. This has been the case in some EU countries before the reform (so-called “C1 beets”) and it is also current practice (see, for example [[http:// | ||
+ | |||
+ | Considering a kinked demand curve and in addition yield uncertainty renders the standard profit maximisation hypothesis inappropriate for the sugar sector (at least). The CAPRI system therefore applies an expected profit maximisation framework that takes care for yield uncertainty (see Adenäuer 2005). The idea behind this is that observed C sugar productions in the past are unlikely to be an outcome of competitiveness at C beet prices rather than being the result of farmers’ aspirations to fulfil their quota rights even in case of a bad harvest. This approach essentially assumes that the “behavioural quotas” of farmers may exceed the “legal quotas” (derived from the sugar CMO) by some percentage. This percentage reflects in part the pricing behaviour of the regional sugar industry, but it may also depend on farmers expectations on the consequences of an incomplete quota fill. These aspects may be captured with the following specification of expected sugar beet revenues that substitute for the expression \(NETTRD_{r, | ||
+ | |||
+ | |||
+ | \begin{align} | ||
+ | \begin{split} | ||
+ | SegbREV_r & = p^A NETTRD_{r, | ||
+ | & - \left( p^A-p^B\right) \left [ | ||
+ | \begin{matrix} | ||
+ | (1-CDFSugb(q^A))(NETTRD_{r, | ||
+ | + (\sigma^S)^2 PDFSugb(q^A) | ||
+ | \end{matrix} | ||
+ | \right] \\ | ||
+ | & -\left(p^B-p^C\right) \left [ | ||
+ | \begin{matrix} | ||
+ | (1-CDFSugb(q^A+B))(NETTRD_{r, | ||
+ | + (\sigma^S)^2 PDFSugb(q^A+B) | ||
+ | \end{matrix} | ||
+ | \right] | ||
+ | |||
+ | \end{split} | ||
+ | \end{align} | ||
+ | |||
+ | Where \(PDFSugb_r\) and \(CDFSugb_r\) are the probability res. cumulated density functions of the NETTRD variable with the standard deviation \(\sigma^S\). \(\sigma^S\) is defined as \(NETTRD_{r, | ||
+ | |||
+ | |||
+ | \begin{align} | ||
+ | \begin{split} | ||
+ | p^A &= legalqout^A \cdot scalefac \\ | ||
+ | & = legalqout^A \cdot \left(\frac{NETTRD_{SUGB}^{cal}}{legalqout^A} \right )^{0.8} | ||
+ | \end{split} | ||
+ | \end{align} | ||
+ | |||
+ | The scaling factor to map from the legal quota legalquotA (as the B quota has been eliminated in the sugar reform, it holds that \(q^A = q^{A+B}) \)to the behavioural quota qA depends on the projected sugar beet sales quantity in the calibration point \(NETTRD_{SUGB}^{cal} : For a country with a high over quota production (say 40%) we would obtain a scaling factor of 1.31, such that this producer will behave like a moderate C-sugar producer: responsive to both the C-beet prices as well as to the quota beet price (and the legal quotas). Without this scaling factor, producers with significant over quota p | ||
+ | |||
+ | ===Update note=== | ||
+ | |||
+ | A number of recent developments are not covered in the previous exposition of supply model equations | ||
+ | |||
+ | -A series of projects have added a distinction of rainfed and irrigated varieties of most crop activities which is the core of the so-called “CAPRI-water” version of the system((A more complete presentation is given in [[https:// | ||
+ | -Several projects have added endogenous GHG mitigation options((These are most completely included in the “trunk” version of the CAPRI system. For details, see, for example, [[http:// | ||
+ | -Several new equations serve to explicitly represent environmental constraints deriving from the Nitrates Directive and the NEC directive((These are most completely included in the “trunk” version of the CAPRI system but developments are still ongoing.)). | ||
+ | -A complete area balance monitoring the land use changes according to the six UNFCCC land use types (cropland, grassland, forest land, wetland, settlements, | ||
+ | |||
+ | ====Calibration of the regional programming models==== | ||
+ | |||
+ | Since the very first CAPRI version, ideas based on Positive Mathematical Programming were used to achieve perfect calibration to observed behaviour – namely regional statistics on cropping pattern, herds and yield – and data base results as the input or feed distribution. The basic idea is to interpret the ‘observed’ situation as a profit maximising choice of the agent, assuming that all constraints and coefficients are correctly specified with the exemption of costs or revenues not included in the model. Any difference between the marginal revenues and the marginal costs found at the base year situation is then mapped into a non-linear cost function, so that marginal revenues and costs are equal for all activities. In order to find the difference between marginal costs and revenues in the model without the non-linear cost function, calibration bounds around the choice variables are introduced. | ||
+ | |||
+ | The reader is now reminded that marginal costs in a programming model without non-linear terms comprise the accounting cost found in the objective and opportunity costs linked to binding resources. The opportunity costs in turn are a function of the accounting costs found in the objective. It is therefore not astonishing that a model where marginal revenues are not equal to marginal revenues at observed activity levels will most probably not produce reliable estimates of opportunity costs. The CAPRI team responded to that problem by defining exogenously the opportunity costs of two major restrictions: | ||
+ | |||
+ | ====Estimating the supply response of the regional programming models==== | ||
+ | |||
+ | The development, | ||
+ | |||
+ | The two possible competitors are standard duality based approaches with a following calibration step or estimates based directly on the Kuhn-Tucker conditions of the programming models. Both may or may not require a priori information to overcome missing degrees of freedom or reduce second or higher moments of estimated parameters. The duality based system estimation approach has the advantage to be well established. Less data is required for the estimation, typically prices and premiums and production quantities. That may be seen as advantage to reduce the amount of more or less constructed information entering the estimation, as input coefficients. However, the calibration process is cumbersome, and the resulting elasticities in simulation experiments will differ from the results of the econometric analysis. | ||
+ | |||
+ | The second approach – estimating parameters using the Kuhn-Tucker-conditions of the model – leads clearly to consistency between the estimation and simulation framework. However, for a model with as many choice variables as CAPRI that straightforward approach may require modifications as well, e.g. by defining the opportunity costs from the feed requirements exogenously. | ||
+ | |||
+ | The dissertation work of Torbjoern Jansson (Jansson 2007) focussed on estimating the CAPRI supply side parameters. The results have been incorporated in the current version. The milk study (2007/08) contributed additional empirical evidence on marginal costs related to milk production (Kempen et al. 2011) | ||
+ | |||
+ | ====Price depending crop yields and input coefficients ==== | ||
+ | |||
+ | Let Y denote yields and j production activities. Yield react via iso-elastic functions to changes in output prices | ||
+ | |||
+ | \begin{equation} | ||
+ | log(Y_j)=\alpha_j+\epsilon_j \, log(p_o) | ||
+ | \end{equation} | ||
+ | |||
+ | The current implementation features yield elasticities for cereals chosen as 0.3, and for oilseeds and potatoes chosen as 0.2. These estimates might be somewhat conservative when compared e.g. with Keeney and Hertel (2008). However, in CAPRI they relate to small scale regional units and single crops, and to European conditions which might be characterized by a combination of higher incentive for extensive management practises and dominance of rainfed agriculture where water might be a yield limiting factor. | ||
+ | |||
+ | Currently, the code is set up as to only capture the effect of output prices. However, in order to spare calculation of the constant terms α, the actual code implemented in ‘// | ||
+ | |||
+ | \begin{equation} | ||
+ | Y_{j, | ||
+ | \end{equation} | ||
+ | |||
+ | ==== LULUCF in the supply model of CAPRI ==== | ||
+ | |||
+ | === Introduction === | ||
+ | |||
+ | This technical paper explains how the most aggregate level of the CAPRI area allocation in the context of the supply models has been re-specified in the TRUSTEE ((https:// | ||
+ | )) and SUPREMA ((https:// | ||
+ | |||
+ | During the subsequent period, CAPRI was increasingly adapted to analyses of greenhouse gas (GHG) emission studies. Examples include CAPRI-ECC, GGELS, ECAMPA-X, AgCLim50-X, (European Commission, Joint Research Centre), ClipByFood (Swedish Energy Board), SUPREMA (H2020). This vein of research is very likely to gain in importance in the future. | ||
+ | |||
+ | In order to improve land related climate gas modelling within CAPRI, it was deemed appropriate to (1) extend the land use modelled to //all// available land in the EU (i.e. not only agriculture), | ||
+ | )), but as always, an operational version emerged only after integrating efforts by researchers in several projects working at various institutions. Within the SUPREMA project another important change in the depiction of land use change was made: the Markov chain approach was replaced by prespecifying the total land transitions as average transitions per year times the projection. This paper focusses on the theory applied while data and technical implementation are only briefly covered. | ||
+ | |||
+ | |||
+ | === A simple theory of land supply === | ||
+ | |||
+ | Recall the dual methodological changes attempted in this paper: | ||
+ | |||
+ | - Extend land use modelling to the entire land area, and | ||
+ | - Explicitly model transitions between each pair of land uses | ||
+ | |||
+ | In order to keep things as simple as possible, we opted for a theory where the decision of how much land to allocate to each use is independent of the explicit transitions between classes. This separation of decisions is simplifying the theoretical derivations, | ||
+ | |||
+ | The land supply and transformation model developed here is a bilevel optimization model. At the higher level (sometimes termed the //outer problem//), the land owner decides how much land to allocate to each aggregate land use based on the rents earned in each use and a set of parameters capturing the costs required in order to ensure that the land is available to the intended use. At the lower level (sometimes termed the //inner problem//), the transitions between land classes are modelled, with the condition that the total land needs of the outer problem are satisfied. The inner problem is modelled as a stochastic process involving no explicit economic model. | ||
+ | |||
+ | For the outer problem, i.e. the land owner’s problem, we propose a quadratic objective function that maximizes the sum of land rents minus a dual cost function. The parameters of the dual cost function were specified in two steps: | ||
+ | |||
+ | - A matrix of land supply elasticities was estimated (by TRUSTEE partner Jean Saveur Ay, CESEAR, Dijon (JSA). This estimation might be updated in future work or replaced with other sources for elasticities. | ||
+ | - The parameters of the dual cost function are specified so that the supply behaviour replicates the estimated elasticities as closely as possible while exactly replicating observed/ | ||
+ | |||
+ | The model is somewhat complicated by the fact that land use classes in CAPRI are defined somewhat differently compared to the UNFCCC accounting and also in the land transition data set. Therefore, some of the land classes used in the land transitions are different from the ones used in the land supply model. In particular, “Other land”, “Wetlands” and “Pasture” are differently defined. To reconcile the differences, | ||
+ | |||
+ | === Inner model – transitions === | ||
+ | |||
+ | A vector of supply of land of various types could result from a wide range of different transitions. The inner model determines the matrix of land transitions that is “most likely”. The concept of “most likely” is formalized by assuming a joint density function for the land transitions, | ||
+ | |||
+ | == Gamma density == | ||
+ | |||
+ | Since each transition is non-negative, | ||
+ | |||
+ | {{: | ||
+ | |||
+ | Figure 1: Gamma density graph for mode=1 and various standard deviations. “acc”=" | ||
+ | |||
+ | Let $i$ denote land use classes in CAPRI definition, whereas //l// and //k// are land uses in UNFCCC classification. Let $\text{LU}_{k}$ be total land use after transitions and $\text{LU}_{l}^{\text{initial}}$ be land use before transitions. Furthermore, | ||
+ | |||
+ | $${\max_{T_{\text{lk}}}{\log{\prod_{\text{lk}}^{}{f\left( T_{\text{lk}}|\alpha_{\text{lk}}, | ||
+ | |||
+ | $$\Rightarrow \max_{T_{\text{lk}}}\sum_{\text{lk}}^{}\left\lbrack \left( \alpha_{\text{lk}} - 1 \right)\log T_{\text{lk}} - \beta_{\text{lk}}T_{\text{lk}} \right\rbrack$$ | ||
+ | |||
+ | subject to | ||
+ | |||
+ | $$\text{LU}_{k} - \sum_{l}^{}T_{\text{lk}} = 0 \; \left\lbrack \tau_{k} \right\rbrack$$ | ||
+ | |||
+ | $$\text{LU}_{l}^{\text{initial}} - \sum_{k}^{}T_{\text{lk}} = 0\; | ||
+ | |||
+ | $$\text{LU}_{k} - \sum_{i}^{}{\text{shar}e_{\text{ki}}\text{LEV}L_{i}} = 0$$ | ||
+ | |||
+ | The last equation is needed to convert land use in UNFCCC classification to land use in CAPRI classification, | ||
+ | |||
+ | $$\ \left( \alpha_{\text{lk}} - 1 \right)T_{\text{lk}}^{- 1} - \beta_{\text{lk}} + \tau_{k}^{} + \tau_{l}^{\text{initial}} = 0$$ | ||
+ | |||
+ | The parameters $\alpha$ and $\beta$ of the gamma density function were computed by assuming that (i) the observed transitions are the mode of the density, and (ii) the standard deviation equals the mode. Then the parameters are obtained by solving the following quadratic system: | ||
+ | |||
+ | $$\text{mode} = \frac{\alpha - 1}{\beta}$$ | ||
+ | |||
+ | $$\text{variance} = \frac{\alpha}{\beta^{2}}$$ | ||
+ | |||
+ | == Annual transitions via Marcov chain in basic model == | ||
+ | |||
+ | The implementation in CAPRI differs from the above general framework in that it explicitly identifies the //annual// transitions in year t $T_{\text{lk}}^{t}$ from the initial $\text{LU}_{l}^{\text{initial}}$ land use to the final land use $\text{LU}_{k}$. This is necessary to identify the annual carbon effects occurring only in the final year in order to add them to the current GHG emissions, say from mineral fertiliser application in the final simulation year. If the initial year is the base year = 2008 and projection is for 2030, then the carbon effects related to the change from the 2008 $\text{LU}_{l}^{\text{initial}}$ to the final land use $\text{LU}_{k}$ (=$T_{\text{lk}}$in the above notation, without time index) refer to a period of 22 years that cannot reasonably be aggregated with the “running” non-CO2 effects from the final year 2030. Furthermore the historical time series used to determine the mode of the gamma density for the transitions also refer to annual transitions. | ||
+ | |||
+ | Initially the problem to link total to annual transitions has been solved by assuming a linear time path from the initial to the final period, but this was criticised as being an inconsistent time path (by FW). Ultimately the time path has been computed therefore in the supply model in line with a static Markov chain with constant probabilities $P_{\text{lk}}$ such that both land use $\text{LU}_{l}^{t}$ as well as transitions $T_{\text{lk}}^{t}$ in absolute ha require a time index (e_luOverTime in supply_model.gms). | ||
+ | |||
+ | $$\text{LU}_{k}^{t} - \sum_{l}^{}{P_{\text{lk}}\text{LU}_{l}^{t - 1}} = 0\ ,\ t = \{ 1,\ldots s\}$$ | ||
+ | |||
+ | Where $\text{LU}_{k}^{s}$ is the final land use in the simulation year s and $\text{LU}_{k}^{0} = \text{LU}_{k}^{\text{iniital}}$ is the initial land use. The transitions in ha in any year may be recovered from previous years land use and the annual (and constant) transition probabilities (e_LUCfromMatrix in supply_model.gms). | ||
+ | |||
+ | $$T_{\text{lk}}^{t} = P_{\text{lk}}*\text{LU}_{l}^{t - 1}$$ | ||
+ | |||
+ | The absolute transitions may enter the carbon accounting (ignored here) and if we substitute the last period’s transitions we are back to the condition for consistent land balancing in the final period from above: | ||
+ | |||
+ | $$\text{LU}_{k}^{s} = \sum_{l}^{}{P_{\text{lk}}\text{LU}_{l}^{s - 1}} = \sum_{l}^{}T_{\text{lk}}^{s}$$ | ||
+ | |||
+ | When using the transition probabilities in the consistency condition for initial land use we obtain | ||
+ | |||
+ | $$\text{LU}_{l}^{\text{initial}} - \sum_{k}^{}T_{\text{lk}}^{1} = 0$$ | ||
+ | |||
+ | $$\Longleftrightarrow \text{LU}_{l}^{\text{initial}} = \sum_{k}^{}{P_{\text{lk}}^{}\text{LU}}_{l}^{\text{iniital}}$$ | ||
+ | |||
+ | $$\Leftrightarrow 1 = \sum_{k}^{}P_{\text{lk}}$$ | ||
+ | |||
+ | So the simple condition is that probabilities have to add up to one (e_addUpTransMatrix in supply_model.gms). | ||
+ | |||
+ | == Annual transitions if SUPREMA is active == | ||
+ | |||
+ | As the use of the Marcov-chain approach allows the annual transitions to be explicit model variables that could be used to compute annual carbon effects but leads to computational limitations especially in the market model a new approach was developed under SUPREMA (i.e. if %supremaSup% == on) by re-specifying the total land transitions as average transitions per year times the projection horizon and by considering for the remaining class without land use change (on the diagonal of the land transition matrix) only the annual carbon effects per ha, relevant for the case of gains via forest management. | ||
+ | |||
+ | The new accounting in the CAPRI global supply model may be explained as follows, starting from a calculation of the total GHG effects G over horizon h = t-s from total land transitions L< | ||
+ | |||
+ | $$G = Γ*h = \sum_{i, | ||
+ | |||
+ | Where Γ collects the annual GHG effects that correspond to the total GHG effects divided by the time horizon G / h. These annual effects may be calculated as based on average annual transitions and annual effects for the remaining class as follows: | ||
+ | |||
+ | $$Γ= \sum_{i, | ||
+ | + \sum_{i}^{}{ε_{\text{ii}}^{}\text{L}}_{ii}^{} $$ | ||
+ | |||
+ | Where Λ< | ||
+ | |||
+ | Using these average annual transitions for true (off-diagonal) LUC we may compute the final classes as follows: | ||
+ | |||
+ | $$ \text{LU}_{k, | ||
+ | |||
+ | While adding up of shares (or probabilities) of LUC from class I to k over all receiving classes k continues to hold as stated above. It should be highlighted that the land use accounting implemented under SUPREMA avoids the need to explicitly trace the annual transitions in the form of a Markov chain and thereby economised on equations and variables. | ||
+ | |||
+ | |||
+ | |||
+ | === Outer model – land supply === | ||
+ | The outer problem is defined as a maximization of the sum of land rents minus a quadratic cost term, subject to the first order optimality conditions of the inner problem: | ||
+ | |||
+ | $$\max{\sum_{i}^{}{\text{LEV}L_{i}r_{i}} - \sum_{i}^{}{\text{LEV}L_{i}c_{i}} - \frac{1}{2}\sum_{\text{ij}}^{}{\text{LEV}L_{i}D_{\text{ij}}\text{LEV}L_{j}}}$$ | ||
+ | |||
+ | subject to, | ||
+ | |||
+ | $$\text{LU}_{k} - \sum_{i}^{}{\text{shar}e_{\text{ki}}\text{LEV}L_{i}} = 0$$ | ||
+ | |||
+ | $$\text{LU}_{k} - \sum_{l}^{}T_{\text{lk}} = 0\; | ||
+ | |||
+ | $$\text{LU}_{l}^{\text{initial}} - \sum_{k}^{}T_{\text{lk}} = 0\; | ||
+ | |||
+ | $$\ \left( \alpha_{\text{lk}} - 1 \right)T_{\text{lk}}^{- 1} - \beta_{\text{lk}} + \tau_{k}^{} + \tau_{l}^{\text{initial}} = 0$$ | ||
+ | |||
+ | The parameters of the inner model **α** and **β// | ||
+ | |||
+ | There are a few methodological and numerical challenges to overcome. In particular, we need to (i) analytically derive $\mathbf{\eta}\left( \mathbf{c}, | ||
+ | |||
+ | $$\sum_{i}^{}{\text{LEV}L_{i}} - \sum_{l}^{}{LU_{l}^{\text{initial}}} = 0$$ | ||
+ | |||
+ | Note that the second sum is a constant. This simplification is based on the observation that the land transitions don’t appear in the objective function of the outer problem, so that all solutions to the inner problems are equivalent from the perspective of the outer problem, and that any land use vector that preserves the initial land endowment is a feasible solution to the inner problem. | ||
+ | |||
+ | Next, we formulate the first order condition (FOC) of the modified outer problem to obtain land use as an implicit function of the parameters, $F\left( LEVL, | ||
+ | |||
+ | The first order conditions, and the implicit function, become | ||
+ | |||
+ | $$F\left( LEVL, | ||
+ | \frac{\partial\mathcal{L}}{\partial LEVL_{i}} = & r_{i} - c_{i} - \sum_{j}^{}{D_{\text{ij}}\text{LEV}L_{j}} - \lambda & = 0 \\ | ||
+ | \frac{\partial\mathcal{L}}{\partial\lambda} = & \sum_{i}^{}{\text{LEV}L_{i}} - \sum_{l}^{}{LU_{l}^{\text{initial}}} & = 0 \\ | ||
+ | \end{bmatrix}$$ | ||
+ | |||
+ | In order to apply the implicit function theorem((Recall that the implicit function theorem states that if F(x,p) = 0, then dx/dp = -[dF/ | ||
+ | )) we need to differentiate the FOC once w.r.t. the variables $\text{LEV}L_{i}$ and $\lambda$ and once with respect to the parameter of interest, $r_{j}$, invert the former and take the negative of the matrix product. If (currently) irrelevant parameter are omitted, the following matrix of $(N + 1) \times (N + 1)$ is obtained (the “+1” is the uninteresting derivative of total land rent $\lambda$ with respect to individual land class rent $r_{i}$) | ||
+ | |||
+ | $$\left\lbrack \frac{\partial LEVL}{\partial r} \right\rbrack = - \left\lbrack D_{LEVL, | ||
+ | |||
+ | $$\begin{bmatrix} | ||
+ | \frac{\partial LEVL}{\partial r} \\ | ||
+ | \frac{\partial\lambda}{\partial r} \\ | ||
+ | \end{bmatrix} = - \begin{bmatrix} | ||
+ | \frac{\partial F}{\partial LEVL} & \frac{\partial F}{\partial\lambda} \\ | ||
+ | \end{bmatrix}\left\lbrack \frac{\partial F}{\partial r} \right\rbrack$$ | ||
+ | |||
+ | Carrying out the differentiation specifically for land rent // | ||
+ | |||
+ | $$\begin{bmatrix} | ||
+ | \frac{\partial LEVL_{i}}{\partial r_{j}} \\ | ||
+ | \frac{\partial\lambda}{\partial r_{j}} \\ | ||
+ | \end{bmatrix} = - \begin{bmatrix} | ||
+ | \left\lbrack {- D}_{\text{ij}} \right\rbrack & - 1 \\ | ||
+ | - 1' & 0 \\ | ||
+ | \end{bmatrix}^{- 1}\begin{bmatrix} | ||
+ | I \\ | ||
+ | 0 \\ | ||
+ | \end{bmatrix}$$ | ||
+ | |||
+ | Discarding the last row of the resulting $(N + 1) \times N$ matrix finally lets us compute the elasticity as | ||
+ | |||
+ | $$\left\lbrack \eta_{\text{ij}} \right\rbrack = \left\lbrack \frac{\partial LEVL_{i}}{\partial r_{j}} \right\rbrack\left\lbrack \frac{r_{j}}{\text{LEV}L_{i}} \right\rbrack$$ | ||
+ | |||
+ | In the estimation, we assumed that the prior elasticity matrix is the mode of a density where each entry were independently distributed. Furthermore, | ||
+ | |||
+ | $$\max_{\eta, | ||
+ | |||
+ | subject to | ||
+ | |||
+ | $$\left\lbrack \frac{\partial LEVL_{i}}{\partial r_{j}} \right\rbrack = - \begin{bmatrix} | ||
+ | \left\lbrack {- D}_{\text{ij}} \right\rbrack & - 1 \\ | ||
+ | - 1' & 0 \\ | ||
+ | \end{bmatrix}^{- 1}\begin{bmatrix} | ||
+ | I \\ | ||
+ | 0 \\ | ||
+ | \end{bmatrix}$$ | ||
+ | |||
+ | $$\left\lbrack \eta_{\text{ij}} \right\rbrack = \left\lbrack \frac{\partial LEVL_{i}}{\partial r_{j}} \right\rbrack\left\lbrack \frac{r_{j}}{\text{LEV}L_{i}} \right\rbrack$$ | ||
+ | |||
+ | $$\begin{matrix} | ||
+ | & r_{i} - c_{i} - \sum_{j}^{}{D_{\text{ij}}\text{LEV}L_{j}} - \lambda & = 0 \\ | ||
+ | & \sum_{i}^{}{\text{LEV}L_{i}} - \sum_{l}^{}{LU_{l}^{\text{initial}}} & = 0 \\ | ||
+ | \end{matrix}$$ | ||
+ | |||
+ | and the curvature constraint using a stricter variant of the Cholesky factorization | ||
+ | |||
+ | $$D_{\text{ij}}\left( 1 - \delta I_{\text{ij}} \right) = \sum_{k}^{}{U_{\text{ki}}U_{\text{kj}}}$$ | ||
+ | |||
+ | where $\delta$ is a small positive number and $I_{\text{ij}}$ entries of the identity matrix such that the factor $(1 - \delta I_{\text{ij}})$ shrinks the diagonal of the D-matrix, ensuring //strict// positive definiteness instead of // | ||
+ | |||
+ | ==Prior elasticities and area mappings== | ||
+ | |||
+ | The empirical evidence obtained in the TRUSTEE project applied to prior elasticities for land categories based on Corine Land Cover (CLC) data. These categories are also covered in the CAPRI database based on various sources (see the database section in the CAPRI documentation): | ||
+ | |||
+ | The introduction has mentioned already three systems of area categories that need to be distinguished. The first one is the set of area aggregates with good coverage in statistics that has been investigated recently by JS Ay (2016), in the following “JSA”: | ||
+ | |||
+ | $$\text{LEVL} = \left\{ \text{ARAC}, | ||
+ | |||
+ | Where | ||
+ | |||
+ | ARAC = arable crops | ||
+ | |||
+ | FRUN = perennial crops | ||
+ | |||
+ | GRAS = permanent grassland | ||
+ | |||
+ | FORE = forest | ||
+ | |||
+ | ARTIF = artificial surfaces (settlements, | ||
+ | |||
+ | OLND = other land | ||
+ | |||
+ | The above categories are matching reasonably well with the definitions in JSA. A mismatch exists in the classification of paddy (part of ARAC in CAPRI but in the perennial group in JSA) and terrestrial wetlands (part of OLND in CAPRI and a separate category in JSA). Inland waters are considered exogenous in CAPRI and hence not included in the above set LEVL. | ||
+ | |||
+ | For carbon accounting we need to identify the six LU classes from IPCC recommendations and official UNFCCC reporting: | ||
+ | |||
+ | $$LU = \left\{ \text{CROP}, | ||
+ | |||
+ | which is typically indexed below with “l” or “k” ∈ LU and where | ||
+ | |||
+ | CROP = crop land (= sum of arable crops and perennial crops) | ||
+ | |||
+ | GRSLND = grassland in IPCC definition (includes some shrub land and other “nature land”, hence GRSLND> | ||
+ | |||
+ | WETLND = wetland (includes inland waters but also terrestrial wetlands) | ||
+ | |||
+ | RESLND = residual land is that part of OLND not allocated to grassland or wetland, hence RESLND< | ||
+ | |||
+ | FORE = forest | ||
+ | |||
+ | ARTIF = artificial surfaces | ||
+ | |||
+ | In the CAPRI database, in particular for its technical base year, we have estimated an allocation of other land OLND into its components attributable to the UNFCCC classes GRSLND, | ||
+ | |||
+ | $$\text{OLND}^{0} = {\text{OLND}G}^{0} + {\text{OLND}W}^{0} + {\text{OLND}R}^{0}$$ | ||
+ | |||
+ | Lacking better options to make the link between sets LEVL (activity level aggregates) and LU (UNFCCC classes, technically in CAPRI code: set “LUclass”) we will assume that these shares are fixed and may estimate the “mixed” LU areas from activity level aggregates as follows | ||
+ | |||
+ | ^// | ||
+ | |WETLND | ||
+ | |RESLND | ||
+ | |||
+ | which means that the mapping from set LEVL to set LU only uses some fixed shares of LEVL areas that are mapped to a certain LU: | ||
+ | |||
+ | $$LU_k=\sum_i{\text{share}_{\text{i, | ||
+ | |||
+ | where 0 ≤ // | ||
+ | |||
+ | ===Technical implementation=== | ||
+ | |||
+ | The key equations corresponding to the approach explained above are collected in file supply_model.gms or the included files supply/ | ||
+ | |||
+ | // | ||
+ | |||
+ | At this point, it should also be explained that rents for non-agricultural land types were entirely based on assumptions (a certain ratio to agricultural rents). As there were no plans to run scenarios with modified non-agricultural rents, these land rents //r// used in calibration for those land types were subtracted from the “c-paramter”, | ||
+ | |||
+ | Furthermore, | ||
+ | |||
+ | More detailed explanations on the technical implementation are covered elsewhere, for example in the “Training material” included in the EcAMPA-4 deliverable D5. | ||
+ | |||
+ | Concerning the improvements made under SUPREMA from a technical perspective, | ||
+ | |||
+ | === Emission Equations === | ||
+ | |||
+ | Under EcAMPA 3 and partly in earlier projects (inter alia EcAMPA 2) new modelling outputs have been developed for indicators without matching reporting infrastructure helping users to organise the additional information. This applied for example to | ||
+ | |||
+ | 1) Additional CAPRI results on land use results related to the complete area coverage, mappings to UNFCCC area categories and their transitions; | ||
+ | |||
+ | 2) The carbon effects linked to these land transitions. | ||
+ | |||
+ | Furthermore, | ||
+ | |||
+ | The scenarios including the emission equations are only run if %ghgabatement% == on, otherwise emissions are only calculated and not simulated. | ||
+ | |||
+ | The following emission equations have been implemented: | ||
+ | |||
+ | ^**Code** | ||
+ | |GWPA |Agricultural emissions | ||
+ | |CH4ENT | ||
+ | |CH4MAN | ||
+ | |CH4RIC | ||
+ | |N2OMAN | ||
+ | |N2OAPP | ||
+ | |N2OGRA | ||
+ | |N2OSYN | ||
+ | |N2OCRO | ||
+ | |N2OAMM | ||
+ | |N2OLEA | ||
+ | |N2OHIS | ||
+ | |GLUC |Emissions related to indirect land use changes | ||
+ | |CO2BIO | ||
+ | |CO2SOI | ||
+ | |CO2HIS\\ \\ CH4HIS|Carbon dioxide emissions from the cultivation of histosols\\ \\ Methane emissions from cultivation of histosols| | ||
+ | |CO2LIM\\ \\ CO2BUR|Carbon dioxide emissions from limestone and dolomit\\ \\ Carbon dioxide emissions from burning | ||
+ | |CH4BUR | ||
+ | |N2OBUR | ||
+ | |N2OSOI | ||
+ | |GPRD |Emissions related to the production of non-agricultural inputs to agriculture | ||
+ | |N2OPRD | ||
+ | |O2PRD | ||
module_for_agricultural_supply_at_regional_level.1583923590.txt.gz · Last modified: 2022/11/07 10:23 (external edit)