The Optimal Selective Logging Regime and the Faustmann Formula

This study analyzes the optimal selective logging regime of a size-distributed forest where individual trees compete for scarce resources such as space, light, and nutrients. The decision problem of the forest manager is formulated as a distributed optimal control problem. The interpretation of the ﬁrst-order conditions allows a generalization of the Faustmann formula. In an empirical part, this article numerically determines the optimal management regime of a size-structured forest and shows that the optimal selective logging regime is associated with a normal forest under a wider variety of situations than stated in the previous literature


Introduction
consequently the size of a forest, is usually measured by the diameter at breast height, i. e., the diameter of 118 the trunk at a height of 1.30 m above the ground. We denote the diameter of a tree by l ∈ Ω, Ω ≡ [l 0 , l m ), 119 where l 0 and l m indicate the biological minimum and maximum size of a tree. The exogenous variable l 120 together with calendar time t form the domain of the control and state variables. We assume that a diameter-121 distributed forest can be fully characterized by the number of trees and by the distribution of the diameter 122 of the trees. In other words, space and the local environmental conditions of the trees are not taken into 123 account. Given that the diameter value of a tree lies in the interval [l 0 , l m ), and that the number of trees is 124 large by assumption, the distribution of the trees can be represented by a density function, denoted by x(t, l), 125 which indicates the population density with respect to the structuring variable, l, at time t. Therefore, the 126 number of trees in the forest at time t is given by 127 X(t) = lm l 0 x(t, l) dl.
(1) the stand basal area to take account of the competition between individuals, 4 that is, 139 E(t) = lm l 0 π 4 l 2 x(t, l) d l.
(2) 140 Hence, the change in the diameter over time of a single tree is given by The instantaneous death rate is denoted by δ(E(t), l). It describes the rate at which the probability of 143 survival of a tree, given the environmental characteristics E(t), decreases with time.

144
The reproduction of the forest can be modeled as internal reproduction or planting. In the former case we 145 would obtain a boundary condition for the partial integrodifferential equation that reflects the reproduction 146 process. Since seed production by individual trees is very high (Karlsson andÖrlander, 2000), it is space, 147 light, and nutrients that are the limiting factors for the upgrowth of young trees, and not the reproduction 148 process itself. For this reason we can assume that the number of seeds that turn into seedlings is sufficiently 149 large. This allows the forest manager to choose the number of trees with a diameter of l 0 by removing any 150 additional trees. The number of upgrowing trees chosen is denoted by p(t, l 0 ). In the case of planting, we 151 are dealing with a forest that is completely managed, where young trees with a diameter of l 0 are planted and 152 no biological reproduction takes place. Hence, for both reproduction systems the control variables for the 153 management of the forest are given by u(t, l) and p(t, l 0 ), denoting cutting density at time t with diameter 154 l, and the flux of the trees with diameter l 0 respectively. Young trees are either grown up to diameter l 0 , or 155 planted with diameter l 0 , at time t. Thus, the optimal forest management problem is a distributed optimal 156 control problem where the time dependent control variable u(t, l) is distributed over l (Feichtinger and diameter l 0 of the tree. Although our modelling approach allows consideration of both planting and natural reproduction, we frame our analysis in the context of natural reproduction since employed selective logging 160 regimes frequently rest on natural reproduction. Based on the well known McKendrick equation for age 161 structured populations (McKendrick, 1926), the dynamics of the diameter-distributed forest can be described 162 by the following partial integrodifferential equation discussed by de Roos (1997), or by Metz and Diekmann 163 (1986) 164 ∂x(t, l) ∂t + ∂g(E(t), l) x(t, l) ∂l = −δ(E(t), l) x(t, l) − u(t, l) (4) 165 subject to the boundary condition g(E(t), l 0 ) x(t, l 0 ) = p(t, l 0 ). The two terms on the left-hand side of 166 equation 4 present the change in the tree density over time and diameter; the second term not only considers 167 the diameter but also takes the interdependence between diameter and time ∂(gx) ∂l = g ∂x ∂l + ∂g ∂l x into account, 168 i.e., it presents the temporal change in diameter multiplied by the change in tree density over diameter plus 169 the temporal change in diameter over diameter multiplied by the tree density. 5 Hence, the flux of the tree 170 density with respect to diameter and time has to equal the terms of the right-hand side of equation 4, given 171 by the product of the mortality rate and the tree density, and the density of the logged trees. 173 We assume that the forest is privately owned and managed over a planning horizon of t 1 . Using the defini-174 tions given in the preceding section, the formal decision problem of the forest owner can be stated as:

The distributed optimal control problem
where E(t) is given by equation (2), and r denotes the discount rate. The twice-differentiable function 179 B(x, u)e −rt presents the discounted net benefits of the timber. It depends not only on the number of logged 180 trees but also on the number of standing trees since it incorporates the maintenance cost of the forest. It 181 is assumed that the maintenance cost function is concave, and thus B x < 0, B xx < 0. It is also assumed 182 that the net benefit function is strictly concave in u. The twice-differentiable and strictly convex func-183 tion C(p)e −rt expresses the discounted cost of nursing trees up to diameter l 0 , the differentiable function 184 S t 1 (x)e −rt 1 the discounted value of the standing trees at the final point in time of the planning horizon, and 185 the differentiable function S lm (x)e −rt expresses the discounted value of the standing trees that have reached 186 the maximum diameter l m . 6 The term x 0 (l) denotes the initial diameter distribution of the trees, and the 187 restriction g(E(t), l 0 ) x(t, l 0 ) = p(t, l 0 ) requires that the flux of the change in diameter at diameter l 0 mul-188 tiplied by the tree density coincides with the total flux of the diameter of trees with diameter l 0 . Finally, the 189 control variables must be nonnegative.

190
Let the costate variable related with the dynamics of the forest be denoted by λ(t, l), and the Lagrange mul-191 tiplier related with the restriction g(E(t), l 0 ) x(t, l 0 ) = p(t, l 0 ) by λ 0 (t ∂λ(t, l) ∂t

194
where µ 1 and η are Kuhn-Tucker multipliers related to the non-negativity constraints of the decision vari-195 ables u and p, respectively. For an interior solution the first necessary condition, equation (5) states that 196 along the optimal path the discounted marginal net benefits of timber should equal the shadow price λ (in 197 situ value) of the forest stock for every t and l. In contrast to lumped optimal control, distributed opti-198 mal control requires that this equation holds along the optimal path not only with respect to time, but also 199 with respect to diameter. Thus, the forest manager maximizes his/her benefits not only over time but also sary condition is just a restatement of the law of motion, and therefore, it will not be discussed here. Fi-212 nally, following Sage (1968) the following necessary transversality conditions have to be taken into account.
The first transversality condition, equation (10), requires at every moment that the shadow cost for nursing 215 trees has to equal the shadow price of the stock evaluated at the diameter l 0 . This transversality condition 216 is a result of the link between the distributed and the boundary control formed by their common stock 217 variable. The transversality condition (11) states that the shadow price of the trees has to equal the value 218 of an additional standing "tree" at the end of the planning horizon. Finally, transversality condition (12) 219 yields that the shadow price of the trees has to be equal to the value of an additional standing tree with the 220 maximum diameter.

Proposition 1 (Generalized Faustmann Equation)
231 The change in the in situ value of the stock density 232 ∂λ(t, l) ∂t

234
can be interpreted as a generalization of the Faustmann formula. to be identical along the optimal path. To see that equation (13) reflects the Faustmann formula as a special 244 case let us restate the Faustmann formula given by 246 where T indicates the age when the entire even-aged stand is cut, P w is the market price of the wood, f (T ) 247 the merchantable volume of wood that a stand of age T produces. The parameter c presents transaction 248 (logging, processing, transport) and nursing costs.

249
The left-hand sides of the Faustmann formula and of equation (13) reflect the change in the in situ value of 250 the forest over time. However, since equation (14) has to be interpreted differently from equation (13), we 251 demonstrate in the appendix that in the case of an even-aged stand, the change in the in situ value over time

252
-the equivalent of equation (13)  value from t to t 1 . Specifically, it holds as t 1 tends to infinity and we obtain that the value of λ(t, l) is given , where the second term reflects the opportunity cost of the land for 259 an infinite stand rotation. Hence, the change in the in situ value dλ/dT is given by P w f ′ (T ) − r c e rT and 260 reflects the left-hand side of equation (14), and rλ = r ) the right-hand 261 side of equation (14) taking into account that the term r c e rT cancels out on both sides, i. e., 262 i.e., all derivatives with respect to l and E are zero.

267
In practical terms, it is not possible to determine whether the optimal rotation age for most trees is lower  3.2 Implications of a Steady state on the optimal long-run distribution.

273
The first order conditions of problem (D) do not permit an analytical solution. To take the analysis further 274 we consider the case where the forest has reached the steady state, i. e. ∂λ/∂t = 0 and ∂x/∂t = 0.

275
The optimal steady state problem is a lumped optimal control problem defined over diameter l ∈ [l 0 , l m ]

276
where the density of trees with diameter l is the control variable. This problem can be thought of as a 277 situation where the private owner chooses the optimal diameter distribution of the trees in the steady state 278 resulting from an exogenous shock. Thus, the optimal path of the stock variable x(l) traces out the optimal 279 steady state distribution. The assumption of an steady-state distribution implies that E(t) is constant, i.e.

280
the density effect is constant and E(t) = E. In this way, the integrodifferential equations (8) and (9) are 281 ordinary differential equations and are mathematically tractable.

282
By suppressing t, equation (9) can be written as: In the case where the benefit function is linear in u(l), the system leads to a corner solution, that is, there 285 exists al ∈ (l 0 , l m ) where u(l) = x(l), and x(l) = 0, ∀l >l. Thus, equation (9), for l <l can be written as

289
where we made use of the boundary condition g(E, l 0 )x(l 0 ) = p(l 0 ). Thus, the optimal long-run distribution 290 will be increasing (decreasing) over the diameter l if the absolute value of dg(E, l) d l is greater (lower) than 291 δ(l).

292
For the case of a low mortality rate and fast growing trees, the optimal steady state distribution is increasing 293 in diameter, i.e., the proportion of large trees is relatively high. On the contrary, when the mortality is high 294 and the growth rate is low, the optimal steady state distribution is decreasing in diameter, i.e., the proportion 295 of large trees is relatively low. In practice the necessary conditions (three equations and a system of partial integrodifferential equations) 298 can only be solved analytically under severe restrictions with respect to the specification of the mathematical 299 problem (Muzicant, 1980). Thus, one has to resort to numerical techniques to solve the distributed control 300 problem. Available techniques such as the method of finite differences, the Galerkin method or that of fi-301 nite elements may be appropriate choices (Calvo and Goetz, 2001 to account for the incorporation of decision variables can be found in Goetz et al. (2008). The Escalator

307
Boxcar Train (EBT) is based on converting the structuring variable into a state variable of the system by 308 transforming the partial integrodifferential equation into ordinary differential equations over time. More-309 over, EBT allows the density effect of the biological processes to be considered. In contrast to the other 310 available methods, it can be implemented with standard computer software used to solve mathematical pro-311 gramming problems.

312
The purpose of the empirical analysis is to determine the optimal selective-logging regime of a diameterdistributed forest, i.e., the selective logging regime that maximizes the discounted private net benefits from 314 timber production of a stand of Pinus sylvestris (Scots pine) over a time horizon of 300 years. 7 In this way 315 it is possible to provide guidance for forest practitioners with respect to the optimal logging pattern, and the 316 optimal long-run diameter distribution of the stand. In order to solve the decision problem (D) it needs to be transformed into a problem which can be solved 319 numerically. For this purpose we define i = 0, · · · , n cohorts over the diameter l, i.e., the trees whose 320 diameters fall within the limits l i and l i+1 are grouped in the cohort i. Hence we can define X i (t) as the 321 number of trees, L i (t) as the average diameter, and U i (t) as the number of cut trees within the cohort i.X,

323
The vectorX 0 denotes the initial density of each cohort. As demonstrated in Goetz et al. (2008), the 324 decision problem (D) can be approximated to the decision problem (D') given by subject to the constraints 327 7 In countries at higher latitudes the species Pinus sylvestris is often considered as shade intolerant, and consequently not suitable for a selective logging regime. However, in countries at mid-range latitudes like Spain natural reproduction requires that older trees protect young trees against heat and water stress in the summer (personal communication by C. Gracia, University of Barcelona, Department of Ecology and CREAF, the Centre for Ecological Research and Forestry Applications). This finding is supported by field experiments reported by Clapham et al. (2002) and Sanchez-Gomez et al. (2006).
where p(t, l 0 ) is now written as P (t) to unify the notation. The term environmentẼ(t) is determined by given by: [0, 50] and divide it into 10 initial subintervals of identical length. In this way, each cohort comprises trees that differ in their diameter by a maximum of 5 cm, and can therefore be considered as homogeneous. The 344 initial number of trees in each cohort, X i (0), i = 1, · · · , n, is calculated in such a way that the basal area 345 of the stand is constant (25 m 2 /ha) in order to allow for comparisons between the results of the different 346 optimization outcomes.

347
The function B(X(t),L(t),Ū (t)) accounts for the net benefits of the timber at time t and is defined as: is given by C(P ) = 0.73 P .

364
The value of the parameters of tree volume, tv(L i ), and the marketable part of the tree volume, mv(L i ), 365 are estimated using information provided by a study by Cañellas et al. (2000). The tree volume follows 366 the allometric relation tv(L) = 0.00157387L 1.745087 , and the marketable part of the volume of timber of each tree is an increasing function of the diameter, given by mv(L) = 0.699 + 0.0004311L. The thinning 368 and nursing period, △t, is set equal to 10 years, which is a common practice for a Pinus sylvestris forest 369 (Cañellas et al., 2000).

370
To determine the dynamics of the forest the growth of a diameter-distributed stand of Pinus sylvestris with-371 out thinning was simulated with the bio-physical simulation model GOTILWA (Growth Of Trees Is Limited 372 by Water 8 ). About 100 different simulation runs were conducted by varying the initial diameter distribution.

373
The results of the simulation were used to estimate the function g(E, L i ), which describes the rate of diame-

387
The mathematical optimization problem (D') was programmed in GAMS (General Algebraic Modeling Sys-388 tem) (Brooke et al., 1992). For the numerical solution of this problem the Conopt3 solver, available within 389 GAMS, was employed. For a given initial distribution, the numerical solution of the problem determines 390 for every 10-year period the optimal logging, U i , and planting, P ; the optimal values of the state variables, 391 X i and L i ; and consequently, economic variables, such as the revenue from timber sales and the cutting, Forest managers who want to maximize net timber benefits have to decide on the intensity of cutting, that 398 is, how many trees of diameter L i have to be cut at each moment of time. To calculate the optimal logging 399 regime we assume that the initial diameter distribution of the trees is given by a beta density function with 400 parameters γ = 0.8 and φ = 0.2, corresponding to a young forest distribution. Table 1   To illustrate the optimal evolution of the forest, Figure 1 a)

422
To illustrate how the initial diameter distribution of the trees alters the optimal selective logging regime, 423 problem (D') has also been solved for an old forest distribution (γ = 2, φ = 0.8), a bell shaped distribution 424 (γ = φ = 2), a U-shaped distribution (γ = φ = 0.5), a uniform distribution (γ = φ = 1) and a non-425 structured forest (even-sized, γ = φ = ∞). Figure 2 depicts the optimal evolution of the weighted average 426 (2a) and standard deviation (2b) of the diameter distribution over time for the analyzed initial distributions. 11 427 10 Please note that the average diameter of the bars (cohorts) is not constant over time. This is explained by the fact that the trees of a cohort always stay together and do not move from one cohort to another. However, since the trees grow the average diameter of the cohort increases as the cohort moves along the time axis. 11 Figure 2 shows the development of the average and standard deviation of the diameter distribution of three different initial diameter distributions. The remaining three initial diameter distributions are not depicted because they follow the same pattern, and their graphical representation would obstruct the interpretation of Figure 2.
It shows that the average diameter of the different distributions tends to converge after approximately 200 years, as the amplitude and phase of the cyclical behavior vanishes. Additionally, Figure 2 shows that the 429 standard deviation of the distributions is governed by the same cyclical pattern. When the mean diameter 430 and standard deviation of the initial diameter distribution are close to that of the steady state (young forest 431 with γ = 0.8, φ = 2), the cyclical evolution of these variables is less pronounced, implying that the benefits 432 will be more stable over time. In general, it can be observed that the long-run mean and standard deviation  Moreover, we conducted a sensitivity analysis to determine the effect of a change in the initial basal area on 439 the steady-state distribution. Figure 3 illustrates the optimal evolution of the weighted average diameter (3a) 440 and the standard deviation (3b) of the diameter over time for a young forest, given the initial basal areas of 441 15, 25 and 35 m 2 /ha. One can see from Figure 3 that the long-run mean and standard deviation tend to the 442 same values as in Figure 2 (17 and 9, respectively). This result shows that the steady-state distribution is not 443 only independent from the initial diameter distribution but is also independent from the initial basal area.  We also conducted a sensitivity analysis to evaluate how the optimal management regime of a forest changes 446 as a result of a variation in the discount rate. Thus, we solved problem (D') for a young forest distribution, 447 given discount rates of 3% and 4%. Figure 4 depicts the optimal evolution of the mean diameter (4a) and 448 standard deviation (4b) over time resulting from the optimizations. Figure 4 shows that the discount rate has 449 a significant influence on the optimal selective logging regime. An increase in the discount rate produces, 450 in the long run, a decrease in the average diameter at which the trees are cut, that is, the trees are cut earlier.

451
However, it can be observed in Figure 4b that the initial diameter distributions stabilize in the long run, 452 independently of the chosen discount rate, i.e., the steady state distribution is independent of the discount 453 rate.
454 Figure 3 455 At the end of section 3.2 we stated that the shape of the distribution at the steady state depends on growth 456 and mortality rates. To illustrate this point, the optimal selective logging regime in the case of an initial 457 diameter distribution of trees corresponding to a young forest is calculated for a higher mortality rate of 458 δ = 0.1 compared to the previously chosen mortality rate of 0.01. The histogram of the resulting long-run 459 distribution is depicted in Figure 5. It shows that an increase in the mortality rate causes the optimal steady-460 state distribution to be decreasing in diameter.  Our results show that in the presence of density effects the steady state distribution tends to a normal forest 463 and is independent of the initial distribution of trees. For a comparison of our results with the previous 464 results in the literature we refer primarily to the article by Salo and Tahvonen (2002) since the employed 465 model is still the up-to-date cornerstone for other works. In the results of their work the optimal long-run age 466 distribution is non cyclical only when the length of the discretely measured time period converges towards 467 zero, when the discount rate is zero or when the Faustmann rotation is not unique. Under these conditions 468 a normal forest may result. Otherwise many other outcomes, which do not correspond to a normal forest 469 are possible. Salo and Tahvonen (2002) observe that the analytical framework of a time discrete model is 470 the principal cause for the emergence of logging cycles. The model presented in this article is based on 471 a different analytical framework, which allows relaxing assumptions made by Salo and Tahvonen (2002): 472 the stand is structured with respect to diameter, the economically more relevant variable, and not age; the 473 cleared land does not have to be replanted immediately; the forest manager can vary the number of planted 474 trees per hectare; and it allows partial harvesting of a cohort. Uusivuori and Kuuluvainen (2005) used an 475 analytical framework very similar to Salo and Tahvonen (2002) but allowed also for partial harvesting of 476 the different age classes, leading to long-run distributions being either cyclical or not. According to their 477 findings the noncyclical forest is typically not a normal forest. However, the authors assume in contrast, 478 to this study, and to the article by Salo and Tahvonen (2002), that the price for one cubic feet of timber is 479 independent of the age of the trees. As mentioned above the price of timber per cubic feet usually increases 480 with the diameter (age) of the trees since the wood can be used for the production of more valuable goods.

481
One would expect this assumption to affect the optimal logging regime. For this reason, and due to the fact 482 that trees in the model of Uusivuori and Kuuluvainen (2005) do not only provide timber but also amenity 483 values, it is not possible to compare their results with the results of the cited literature.  The resulting necessary conditions of this problem allow the derivation of an analytical expression which 494 can be interpreted as a generalization of the Faustmann formula. Since the necessary conditions of this 495 problem include a system of partial integrodifferential equations, it usually cannot be solved analytically.

496
Thus, the Escalator Boxcar Train method is proposed to solve the problem numerically. The method allows 497 the partial integrodifferential equation to be transformed into a set of ordinary differential equations and 498 thereby approximate the distributed optimal control problem by a standard optimal control problem. In con-499 trast to the existing literature, the resulting optimization problem can be solved numerically using standard 500 mathematical programming techniques and does not require programming complex numerical algorithms.

501
To determine the optimal selective logging regime of a diameter-distributed and privately owned forest 502 where individual trees compete for scarce resources, an empirical analysis is conducted. It shows that the 503 long-run mean and variance of the diameter distribution for the different types of analyzed forests consid-504 ered tend to a common value, giving rise to a normal forest. Although the diameter distribution of the forests 505 in the steady-state is independent of the initial distribution, the competition between individuals belonging 506 to the same population affects the transition paths to the steady-state distribution, and therefore the opti-507 mal selective-logging regimes in the first periods differ considerably depending on the distribution of the 508 individual characteristics over the entire population.

509
In the case of an even-aged stand decision problem and an infinite planning horizon, problem (D) can be 511 simplified, and is given by 513 subject to the constraints The state and control variables do not depend anymore on the diameter l and since all trees are planted at the 516 same time, the planting costs can be incorporated into the net benefit function of the timber,B(x(t), u(t)).

517
Moreover, the residual value of the stand approaches zero as t approaches infinity, andg(x(t)) denotes As long as I > 0, i.e., we do not cut, µ 2 will be zero and we obtain where v has been replaced by the maximum site value P w f (T ) − c e −rT − 1 . Hence we can conclude that given the 555 optimal value of I, equation (A. 2), which is equivalent to equation (A. 5), provides the Faustmann formula 556 for an even-aged stand.