Integration of an 

Economy under Imperfect 

Competition with a 

Twelve-Cell Ecological Model 
EPA/600/R-06/046 

July 2006 

Integration of an Economy under Imperfect 

Competition with a 

Twelve-Cell Ecological Model

by 

Harland Wm. Whitmore, Jr. 

Department of Economics 

University of Cincinnati 

Cincinnati, OH 45221 

and 

Christopher W. Pawlowski,* Heriberto Cabezas, 

Audrey L. Mayer, and N. Theresa Hoagland 

Sustainable Technology Division 

National Risk Management Research Laboratory 

Cincinnati, OH 45268 

Project Officer 

Heriberto Cabezas 

Sustainable Technology Division 

National Risk Management Research Laboratory 

Cincinnati, OH 45268 

National Risk Management Research Laboratory 

Office of Research and Development 

U.S. Environmental Protection Agency 

Cincinnati, OH 45268 

* Postdoctoral Research Fellow, Oak Ridge Institute for Science & Education (2001-2005) 
Notice 

The U.S. Environmental Protection Agency, through its Office of Research and Development, 

funded and collaborated in the research described here under the Simplified Acquisition Order 

Number 4C-R101-NASA to H.W. Whitmore. It has been subjected to the Agency`s review and 

has been approved for publication as an EPA document. 

Disclaimer 

Mention of trade names or commercial products does not constitute endorsement or 

recommendation for use. 

ii 
Foreword 

The U. S. Environmental Protection Agency (EPA) is charged by Congress with protecting the 

Nation's land, air, and water resources. Under a mandate of national environmental laws, the 

Agency strives to formulate and implement actions leading to a compatible balance between 

human activities and the ability of natural systems to support and nurture life. To meet this 

]mandate, EPA's research program is providing data and technical support for solving 

environmental problems today and building a science knowledge base necessary to manage our 

ecological resources wisely, understand how pollutants affect our health, and prevent or reduce 

environmental risks in the future. 

The National Risk Management Research Laboratory (NRMRL) is the Agency's center for 

investigation of technological and management approaches for preventing and reducing risks 

from pollution that threaten human health and the environment. The focus of the Laboratory's 

research program is on methods and their cost-effectiveness for prevention and control of 

pollution to air, land, water, and subsurface resources; protection of water quality in public water 

systems; remediation of contaminated sites, sediments and ground water; prevention and control 

of indoor air pollution; and restoration of ecosystems. NRMRL collaborates with both public 

and private sector partners to foster technologies that reduce the cost of compliance and to 

anticipate emerging problems. NRMRL's research provides solutions to environmental problems 

by: developing and promoting technologies that protect and improve the environment; advancing 

scientific and engineering information to support regulatory and policy decisions; and providing 

the technical support and information transfer to ensure implementation of environmental 

regulations and strategies at the national, state, and community levels. 

This publication has been produced as part of the Laboratory's strategic long-term research plan. 

It is published and made available by EPA's Office of Research and Development to assist the 

user community and to link researchers with their clients. 

Sally Gutierrez, Director 

National Risk Management Research Laboratory 

iii 
Abstract 

This report documents the scientific research work done to date on developing a generalized 

mathematical model depicting a combined economic-ecological-social system with the goal of 

making it available to the scientific community. The model is preliminary and has not been 

tested or fully explored. The model system described here is intended to represent the first steps 

in combining (in simple fashion) the basic dynamic elements of an ecosystem functioning with a 

human society and an economy in a closed system with a non-limiting supply of energy (the 

model is based on flows of mass between system compartments while the total mass is 

conserved). In this preliminary model, optimizing economic agents (firms and households) 

interact in specific markets and with an ecological system consisting of resource pools and 

several domesticated and wild species. The result is an interdependent system that attempts to 

model macroeconomic variables (based on underlying property rights) and environmental stocks 

and flows. The report contains four chapters. Chapter 1 introduces the project’s objective and 

approach; provides a background on sustainability as a complex system composed of several 

dimensions interacting through time; and provides a brief overview of the integrated 

economic/ecological model and its development. Chapter 2 gives a detailed description of the 

integrated model as it currently exists with equations and explanations. Chapter 3 provides the 

model solution and operational equations. Chapter 4 provides a summary of the report. The 

appendices contain a glossery of terms and the computer code used to implement the model in a 

form suitable for simulating different scenarions. 

A portion of this report was submitted in partial fulfillment of Simplified Acquisition Order 

Number 4C-R101-NASA by H. W. Whitmore, under the sponsorship of the United States 

Environmental Protection Agency. This report covers a period from March 3, 2004 to August 

15, 2005. The work was completed as of February 28, 2006. 

iv 
Table of Contents 

Notice.............................................................................................................................................. ii 

Disclaimer....................................................................................................................................... ii 

Foreword........................................................................................................................................iii 

Abstract.......................................................................................................................................... iv 

Tables............................................................................................................................................. vi 

Figures............................................................................................................................................ vi 

Chapter 1 Introduction ................................................................................................................... 1 

Objective and Approach ............................................................................................................. 1 

Background................................................................................................................................. 2 

Ecological ............................................................................................................................... 4 

Economic................................................................................................................................. 5 

Social/Legal ............................................................................................................................ 6 

Model Overview ......................................................................................................................... 8 

Chapter 2 Integrating the Economic Sector and Ecological Base ............................................... 10 

Economic Decisions of Firms and Households ........................................................................ 11 

Optimal Economic Behavior for P1 Industry ....................................................................... 13 

Optimal Economic Behavior for H1 Industry....................................................................... 17 

Optimal Economic Behavior for IS Industry ........................................................................ 21 

Government Sector................................................................................................................ 24 

GDP ...................................................................................................................................... 24 

Optimal Economic Behavior of the Household Sector ......................................................... 24 

Gross Domestic Product....................................................................................................... 27 

Naturally Occurring Changes in Biological and Physical Resources....................................... 28 

Growth of Non-domesticated Plants, P2 .............................................................................. 28 

Growth of Non-domesticated Herbivores, ............................................................................ 29 

Growth of Non-domesticated Carnivores, C1 ...................................................................... 29 

Growth of Non-domesticated Plants, P3 .............................................................................. 30 

Growth of Non-domesticated Herbivores, H3 ...................................................................... 30 

Growth of Non-domesticated Carnivores, C2 ...................................................................... 30 

Growth in Human Population and Growth in Human Mass ................................................ 31 

Resource Pool ....................................................................................................................... 31 

Inaccessible Resource Pool .................................................................................................. 31 

Chapter 3 Model Solution and Operational Equations ................................................................ 33 

Model Solution.......................................................................................................................... 33 

Operational Equations............................................................................................................... 34 

Chapter 4 Summary ...................................................................................................................... 41 

Appendix A Symbols Used in the Report.................................................................................... 44 

Appendix B SIMULINK Graphical Model ................................................................................. 48 

Appendix C SIMULINK Code .................................................................................................... 49

Appendix D MATLAB Code....................................................................................................... 59 

References..................................................................................................................................... 88 

v 
Tables

Table 3-1 State variables.............................................................................................................. 34 

Table 3-2 Additional model outputs ............................................................................................ 35 

Table 3-3 Parameters and their nominal values........................................................................... 39 

Table 3-4 Corresponding state variable initial conditions. .......................................................... 40 

Figures 

Figure 1-1 Conceptual path of a complex dynamic system........................................................... 3 

Figure 1-2 The integrated model.................................................................................................... 8 

Figure 2-1 The ecomonic model illustrating flows of income and spending .............................. 11 

Figure B-1 Root level SIMULINK model................................................................................... 48 

vi 
Chapter 1 

Introduction 

Objective and Approach 

In many attempts to model economic and ecological systems, the coupling of the 

economic to the ecological in one model is limited, typically interacting through a single 

link or pathway. This method of modeling is inherently inappropriate and perpetuates the 

misconception that human activities can operate outside of the ecological systems on 

which the economic system and human societies depend (O’Neill and Kahn 2000, Rees 

2002). In other cases, economic decision making is integrated into an ecosystem model, 

although the modeling of market mechanisms and the feedbacks between these and the 

ecosystem are limited (e.g., the lake eutrophication and fisheries models of Carpenter et 

al. 1999, Brock and Starrett 2003, Ludwig et al. 2003, and Carpenter and Brock 2004; 

Brock and Xepapadeas 2002; and the Patuxent River watershed models of Voinov et al. 

1999 and Costanza et al. 2002). All of these models were developed to better understand 

and manage human impacts on specific ecological systems; however, they generally have 

very high data requirements for estimating model parameters. 

Several more abstract models have been developed to understand the general behaviors of 

an integrated ecological/economic system, bypassing the need for large amounts of 

ecosystem-specific data. Van den Bergh (1996) offers a basic model that integrates the 

essential features of an ecological system with a fairly simplified production economy. 

His production functions account for materials balance, waste and recycling. In addition, 

he allows for renewable resources to regenerate and for the biosphere to assimilate a 

portion of pollution created during the production process. However, the economic 

portion of his model does not allow for labor services, a pricing mechanism or 

endogenous consumption. The extremely large GUMBO (Global Unified Metamodel of 

the Biosphere) includes an elaborate ecosystem together with production functions and an 

economic welfare function, but no explicit system of market prices (see, for instance, 

Baumans et al. 2002). A wide variety of Computable General Equilibrium (CGE) models 

exist which integrate certain features of the environment into an economic framework; 

they are lucidly reviewed by Conrad (2002). These models do account for the fact that 

emissions and the stock of pollution adversely affect the environment’s ability to provide 

basic services. In addition, they permit environmental quality to influence economic 

efficiency or human welfare. These models contain well developed input-output matrices 

that include material flows across industries. Consumption, production and investment 

are derived endogenously using optimization techniques. Typically however the price of 

labor and the price of capital are exogenously determined. Product prices are determined 

either through perfect competition or by assuming that price is set equal to a constant 

average (and marginal) cost. In the latter case, product demand alone determines the 

level of production. CGE models do not incorporate the dynamics of the regeneration of 

renewable resources, endogenous population growth or the capacity of the biosphere to 

1
assimilate a portion of the pollution generated. 

The preliminary model described here differs in several respects. First, resource limits 

are addressed by explicitly modeling the system as closed to mass; economic decisions 

and biological interactions between species determine how mass is distributed throughout 

the system, and consequentially whether some resources are scarce or some species go 

extinct (i.e., species with zero mass). Second, a legal foundation is incorporated by 

identifying the mass in terms of its property type. Third, an explicit market system of 

decision making is implemented in the form of a price setting model. The model 

introduced here has evolved from much simpler, ecosystem-only models (see Cabezas et 

al. 2003, Fath et al. 2003) and from an ecosystem with economic-like behavior (Cabezas 

et al. 2005). This model also includes an industrial process subsystem in addition to the 

macroeconomic decision-making system and legal foundation. 

The ultimate objective of this modeling effort is to gain general insight as to which legal 

and economic strategies might contribute to or degrade desirable system regimes and 

resilience in an integrated ecological-economic-legal system. It also hoped that some 

understanding is gained about which policies or laws (theoretically) could change 

economic and social behavior to increase the stability of this system. At this point, the 

model is a preliminary one and has not been tested or fully explored. Nor is the model 

calibrated or tied to any real-world systems; rather, it is meant to retain a degree of 

abstractness and generality that may lead to the drawing of general conclusions about the 

interaction between economic markets, ecosystems, and the law. 

Background 

Interest in sustainability has grown exponentially as it has become increasingly obvious 

that the supporting biological systems of the Earth can not indefinitely support current 

rates of human population growth and resource consumption (Millennium Ecosystem 

Assessment Synthesis Reports, 2005). Consider, for example, that according to the 

United States Census Bureau (2005) the human population of the earth grew from 2.5 

billion in 1950 to about 6.4 billion in 2005. Human population growth continues. Also 

consider that from 1970 to 1995, consumption expenditures in 1995 U.S. dollars 

increased from $8.3 to $16.5 trillions in industrialized nations, and from $1.9 to $5.2 

trillions in developing nations (United Nations Development Programme 1998). Lastly 

consider that the human population presently appropriates about 20% of the world net 

terrestrial primary production, leaving a vastly reduced resource pool for all other species 

(Imhoff et al., 2004, Haberl et al. 2004). Here, net primary production is defined as the 

net amount of solar energy used to convert mass to terrestrial plant organic matter by 

photosynthesis. 

Sustainability is fundamentally an effort to create and maintain a regime in which the 

human population and its necessary energy and material consumption can be supported 

indefinitely by the biological system of the Earth. Hence, sustainability is not a goal but 

a path or corridor through time. Figure 1-1 illustrates conceptually the sustainable path 

2
through time of a system through a corridor in a space where the coordinates are 

measurable ecological, industrial, economic and other variables. A sustainability corridor 

is, therefore, defined such that, for example, biodiversity and human population sizes are 

appropriate, the industrial processes perform at high efficiency with minimal 

environmental impacts, and the level of economic activity is adequate to provide 

employment and meet human needs. However, in an integrated system such as this one, 

deviations in any dimension have repercussions elsewhere. For example, inefficient and 

wasteful production causes pollution which damages ecosystems. Therefore, the 

construction of a sustainable corridor requires at least a basic understanding of the 

relationship between production processes, ecosystems, and economies. 

Technological

Economic Dimensions 

Dimensions Catastrophic 

Unsustainable 

Regime 

Self-Correcting Ecological Unsustainable 

Regime Dimensions 

Legal/Social 

Dimensions 

Sustainable 

Regime 

Time 

Figure 1-1 Conceptual path in time for a complex cyclic dynamic system having economic, technological, ecological, legal, and social components. Note the conceptual limits that define sustainable 

regimes shown here as a tunnel, and the three categories of regimes: sustainable, self-correcting unsustainable, and catastrophic unsustainable. 

Hence, as already discussed the definition, assessment, and attainment of sustainability 

are by nature multidisciplinary (Goodland and Daly 1996, Dasgupta et al. 2000, Cabezas 

et al. 2003, McMichael et al. 2003, Cabezas et al. 2005; Figure 1). Popular definitions of 

3
environmental sustainability abound, some of which can be contradictory (World 

Commission on Environment and Development 1987, Gatto 1995, Goodland and Daly 

1996). The overarching concept of sustainability for all disciplines, when applied to 

humanity as a whole or a particular society, highlights the level of activity that can be 

sustained for a given length of time without diminishing the productivity of the system or 

its capacity to recover function following disturbances. While the time frame over which 

this concept is applied can differ markedly between disciplines, all disciplines approach 

sustainability with a reasonably consistent idea of which characteristics of the system are 

desirable, and which are not. From the many possible dimensions of sustainability, three 

which most relevant to the work in this report (ecological, economic, and social/legal) are 

described below. 

Ecological 

Ecological sustainability usually infers that an ecosystem can retain an ability to function 

through environmental changes and disturbances, and over the long term has the 

evolutionary capacity (through genetic and species diversity) to form adaptive ecosystem 

structures and functions. Ecosystems are comprised of collections of species that operate 

in fairly characteristic trophic (feeding) levels, and that engage in a wide variety of 

positive, neutral, and negative interactions. The stability and productivity of ecosystems 

is dictated by the biodiversity of a system, although the role of a specific species or 

population in an ecosystem is not known in all cases with certainty (Tilman 1999, Bond 

and Chase 2002, Hooper et al. 2005). However, due to habitat loss, overharvesting, 

invasion by non-native species, and other anthropogenically induced pressures, many 

species are currently at risk of extinction (Ehrlich 1995). 

As species are lost from systems, the connectivity (and perhaps redundancy) of the 

system declines, and critical functions (such as nutrient recycling, waste treatment or 

pollination) may either no longer be provided (Kearns 1997, Hooper et al. 2005) or are 

provided at reduced capacity. Extinction rates have generally risen proportionately to the 

area of natural habitat impacted by human activity (MacArthur and Wilson 1967, Doak 

and Mills 1994, Pimm and Askins 1995). Extinction (either at the population or species 

level) reduces the evolutionary capacity of ecosystems, and with it their ability to adapt to 

changing environments. Loss of biodiversity in ecosystems would only be sustainable if 

immediate restoration of these systems were possible, but this is unlikely as restoration of 

ecosystems depends heavily on the order of species reintroduced into a system. 

Community assembly rules, resistance to invasion, and the characteristics of each species 

are all critical determinants of restoration efforts, but are at best vaguely understood for 

selected ecosystems (Levin et al. 2001, Sakai et al. 2001, Ferenc et al. 2002). 

Ecosystems are complex, dynamic systems that often display characteristic regimes of 

behavior dictated by their internal dynamics and the disturbances that act on them 

(Scheffer et al. 2001, Mayer and Rietkerk 2004). A dynamic regime, or “alternative 

stable state,” is characterized formally by a particular multidimensional neighborhood or 

set of values over which the system state varies. About each set of such steady states, a 

4
basin of attraction is formed such that absent any changes in external disturbances or 

random variations within the system, the system will remain in that basin of attraction 

and tend toward the steady state. In this analogy, a change in regime is a change to 

another basin of attraction. The size of disturbance that can be tolerated by an ecosystem 

before a change in regime occurs is a measure of its resilience (Holling 1973, 1996, 

Gunderson 2000). Disturbances can range in size, intensity, and frequency, and can 

originate from natural (e.g., lightning strike, fires, floods) or anthropogenic sources (e.g., 

agriculture, deforestation). Although ecosystems may naturally pass through many 

regimes, the functions and services that ecosystems can provide human societies do vary 

under different regimes (Wardle et al. 2000, Portela and Rademacher 2001). In this 

respect, some regimes may be more desirable to humans than other regimes (Carpenter et 

al. 2001). 

Economic 

All economic activity is dependent upon natural resources provided by the environment, 

both as inputs and as sinks for waste mass, e.g., pollutant treatment. Over human time 

scales, these resources can either be renewable (such as fish populations, timber stands, 

or the capacity of the environment to absorb some forms of pollution) or non-renewable 

(such as coal or copper), but renewable resources can be exhausted if harvest rates are too 

high (Slade 1982, Reed 1986, Pauly et al. 2002). Depletion rates of natural resource 

stocks are often dictated by the scarcity of the stock (which influences its price), or the 

economic discount rate used by the industry, instead of the biological limits of 

regeneration of the stock. This discrepancy can lead to extraction and harvest rates for 

renewable resources that are not sustainable. 

Goods and services that ecological systems provide to human economies and societies are 

rarely directly valued in economic markets (Daily 1997). Historically, economists have 

not internalized pollution and other negative impacts to ecosystems in economic analyses, 

but rather have treated these impacts as externalities (Samuelson 1954, Freeman 1984, 

Bird 1987). However, feedback loops between ecological and economic systems can 

significantly alter projected resource availability and economic productivity, and are 

important in the determination of sustainable resource use (Settle et al. 2002). Several 

studies have attempted to make the costs of these services explicit, and estimate that the 

planetary value of all ecosystem services is around $16-54 trillion per year (Costanza et 

al. 1997, although see Heal 2000). Some of these estimates are based on the cost of 

designing and building a technological system that could provide the same services, such 

as reverse osmosis or desalinization technology to produce freshwater in the place of 

wetlands (Postel and Carpenter 1997). Other estimates are based on a variety of hedonic 

pricing, contingent valuation or other more indirect methods (Farber et al. 2002). The 

purpose for valuing ecosystem goods and services is not necessarily to include them in 

regular market transactions, but rather to provide an accounting system that is 

recognizable to economic analysts and can be used to monitor rapid depletion or 

unsustainable use of these goods and services. 

5
In many analyses of economic sustainability, the discount rate chosen and the 

substitutability of environmental capital for other resources and for manmade capital can 

greatly influence the perceived sustainability of an activity. Higher discount rates cause 

firms to extract natural resources faster earlier, saving less for later, and thus leaving a 

reduced (or exhausted) resource base for future generations (Costanza et al. 1997, Barbier 

and Markandya 1998). On the other hand, resources for which there are many 

substitutable sources may be more likely to be used sustainably, as economic agents can 

more easily shift from one resource to another when prices rise (due to increasing 

scarcity). However, as some things have no substitution, such as breathable air or fresh 

water, and economists have long been at odds regarding the degree of substitutability 

between natural and manmade capital (Krutilla 1967). 

Social/Legal 

Since sustainability is ultimately about human well-being, the social aspect in general and 

its legal component in particular must also be considered. Human societies have always 

interacted with and depended upon ecosystems, but the nature of the relationship changed 

dramatically with the onset of agriculture about 10,000 years ago and with the evolving 

concept of property rights. 

The domestication of animals and plants allowed for an exponential increase in human 

population size and density. Food surpluses resulted in job specialization, centralized 

conflict resolution and decision making, collection and redistribution of wealth, 

technology development, and amalgamation of smaller territories into larger ones 

(Diamond 1999). To reflect modern society and the current state of affairs, the plant and 

animal species can be divided into two categories: domesticated and non-domesticated. 

Transfers of mass associated with domesticated species are substantially affected by 

society’s economic system. Transfers of mass associated with non-domesticated species 

are more often affected by the legal and political systems, and/or by biological rules 

alone. Humans have also become adept at appropriating and storing a large portion of the 

nutrient pool by making nutrients unavailable to themselves and the rest of the ecosystem 

for long periods of time (e.g. covering soil with asphalt, water consumed by industry or 

contaminated by pollutants is no longer available for drinking). This is in addition to the 

massive physical infrastructure and energy needed to transport, process, and distribute the 

agricultural products and other resources used to support other societal goals (i.e., 

production of non-food goods and services). 

From a legal standpoint, the ecosystem is less about mass and more about the rights (or 

the lack of rights) associated with the mass. Property is commonly divided into three 

general types: private, state (also referred to as government or public) and commons 

(Heller 2000). Private property has come to be considered a “bundle of rights” that can 

be separated (or “unbundled”) and owned and transferred separately. An example of this 

unbundling, one person (the landlord) may have the ownership rights to a parcel of land, 

while another (the renter) holds the right of possession, and a third (the landlord’s heir) 

has a right to own the parcel at some time in the future. State or government property is 

6
owned or controlled by the sovereign, e.g., a national park. Some legal scholars refer to 

“the commons,” as meaning property which is available to everyone in the world and use 

the term, “common property” to denote a separate category, meaning property available 

to a specific group and excluding others outside that group (Yandle and Morriss 2001). 

Others refer to the former as “open-access commons” and the latter as “closed-access 

commons.” In any case, the difference between the two types of commons is in who has 

the right to use the resource and who has the legal right to exclude others from using it. 

Additional property types have also been identified. Yandle (1999) uses “regulatory 

property” to refer to property which is created and allocated by the government. Air 

pollution permits are an example of this type of property. Heller (1998) describes a type 

of property he calls “anticommons” to account for situations in which property rights are 

so divided among individuals that any one of them can exclude others from effective use 

and the result is underuse of the resource. This is the opposite of open-access commons, 

the well known subject of the Tragedy of the Commons (Hardin 1968), in which no one 

can be excluded and everyone maximizes his/her use of the resource resulting in its 

overuse. 

In many ways, the property rights assigned to the mass in the ecosystem form the basis 

for the economic system because “the process of defining property rights defines wealth 

and its distribution in society” (Yandel and Morriss 2001). Using barbed wire as an 

example, the authors illustrate that it is the transaction costs of defining, defending, and 

devising the property that ultimately determines its fate. Before the invention of barbed 

wire, it was too expensive to enclose and defend large areas of rangeland and so it 

remained open-access. Once the wire was created, it became economically feasible to 

enclose the area (define it) and enforce the exclusion of others (defend it). At that point, 

the rancher had sufficient ownership control to have something worth selling to others 

(devise it). The other option would have been for the rangeland to stay in or revert to 

government ownership with grazing rights sold as “regulatory property.” There are still 

transaction costs involved in this option, but they are borne by taxpayers as well as the 

permit purchaser. 

Transaction costs for open-access commons resources are lower because no rights must 

be defined, defended, or devised, but mass in private hands is more subject to control and 

manipulation for maximum utility. Therefore, the hunter-gatherer has more in common 

with the other mobile species of the foodweb. Both value open-access commons or 

common property because enclosures and exclusive rights pose problems for moving 

across space. The food producer has opposite goals, because enclosures and exclusive 

rights increase wealth--the food producer benefits most by separating its domesticated 

species from others. The modern consumer is aligned with the food producer, since food 

production allows for a much more stable and convenient lifestyle. 

Thus, human society with its legal system of property rights affects ecosystem 

sustainability in several ways. It manipulates foodwebs in favor of domesticated species, 

potentially affecting resources available to non-domesticated species. It appropriates 

portions of the resource pool for physical and social infrastructure. Finally, it raises the 

7
bar on sustainability by a cyclical process of continually increasing its needs through 

population growth made possible by food surpluses brought about by manipulation of the 

food web. 

Model Overview 

The model (Figure 1-2) was developed in stages. The process began by establishing a 

basic ecological model with the total mass distributed among all compartments within a 

foodweb (Fath et al. 2003). The compartments were then differentiated in terms of 

human control over them by identifying species as being either domesticated or nondomesticated, assigning property rights to each compartment, and indicating intentional 

changes to the natural flows of mass, such as fences. Finally, a generic industrial process 

was added to account for resources diverted to human (non-food) consumption and use 

(Cabezas et al. 2005). 

Resource Pool 

Inaccessible Resource Pool 

P1 

HH 

C1 C2

H2 H3

H1 

P2 P3

IS 

Figure 1-2 Integrated model showing flows of mass between compartments and property types. P1, 

H1, and IS are private property. P2 is state-owned property, but H1 has access to it through grazing 

leases issued by the government. C1 is a protected species, thus, H1’s access to it is limited. P3, H2, 

H3, C2 are open-access commons to which no property rights attach. The resource pool in this 

model is open-access commons. Dotted lines indicate mass flows that occur under anthopogenic influence. Gray lines from the inaccessible resource pool to P2 and P3 indicate slow transfers of mass 

as a result of bacterial decay (this is the only natural outlet for mass to escape the inaccessible resource pool). Some mass is also transferred from each compartment to the resource pool representing the death of biological mass; however, to avoid confusing lines, these transfers to RP are not 

shown. 

8
The resulting integrated model comprises twelve compartments, including two resource 

pools (RP and IRP), three primary producers (plants P1, P2, and P3), three herbivores 

(H1, H2, and H3), two carnivores (C1 and C2), an industrial sector (IS) and humans 

(HH). The system flows throughout are specified in terms of mass. The system is closed 

to mass (i.e., mass is conserved) and open to energy. Individual compartments, including 

those composed of private property, observe conservation of mass; that is, any difference 

between input and output in a compartment results in a corresponding change in mass in 

the compartment. Primary producers make available resources from an accessible 

resource pool (RP) to the rest of the food web. Although not shown on Figure 1-2, all 

biological compartments (i.e., all compartment except for IRP and IP) recycle mass back 

to the RP through death. There is a flow of mass from all biological compartments to RP 

proportional to the death rate of the species represented by the compartment. Mass from 

the IRP is recycled very slowly back to P2 and P3 through the action of bacteria; thus, it 

is “inaccessible”to the other compartments. The flow of mass through the biological part 

of the system is determined by a set of Lotka-Volterra type expressions, while the 

resource pools (RP and IRP) simply follow a simple mass balance, and the industrial 

sector (IP) follows a simple flow through with no mass accumulation. 

As discussed, the model structure depicted in Figure 1-2 does not represent any particular 

real ecosystem. Rather it is meant to capture some of the typical features of combined 

ecological-economic-social systems such as: (1) an organization based on trophic levels, 

(2) decreasing number of species with higher tropic levels, (3) species specific 

preferences for food source (e.g., C1 consumes H2 but not H1), (4) the presence of 

humans with industrial production, agricultural production, an economy, and law 

including private property, and (5) the presence of mass that is biologically unavailable as 

result of industrial activity. While the specific structure is arbitrary, it is carefully 

constructed to try to capture as many of these typical features as possible. It is hoped that 

it has enough generality to produce some basic insights into the mechanics of combined 

ecological-economic-social systems. In essence, the model is the analog of a simple 

machine (e.g., a pendulum), which while simple and arbitrary, can still be used to study 

and illustrate basic laws of mechanics that are applicable to more complex machines 

(e.g., a gear box). It is expected that the model, perhaps with some modification, can be 

used to simulate the ecological consequences of different economic and regulatory 

strategies. In this sense this model would be a valuable tool for generically exploring 

sustainable environmental management strategies, but without the risk of experimenting 

with real ecosystems and people. 

9
Chapter 2 

Integrating the Economic Sector and Ecological Base 

This chapter provides the equations and rationale for integrating the economic and 

ecological aspects of the model. It includes the optimal economic behavior of industry, 

government, and households, the natural growth of non-domesticated species and the 

human population, as well as the growth and depletion of the resource pools. 

From an economic perspective the model contains human households (HH), an industrial 

sector (IS), and two private firms: one a producer of plants (P1) and one a producer of 

herbivores (H1). The households are the ultimate owners of the factors of production 

used to produce the goods that are traded in explicit markets. In addition to markets for 

the three goods (P1, H1, IS), there is a labor market. Households must decide between 

devoting time to working for one of the three industries or to leisure. Household income 

comprises labor income and profits generated by the three industries (P1, H1 and IS). In 

this model the stock of physical capital is held constant. A single period planning 

horizon is assumed, therefore savings and investment are ignored. Any dividend income 

is divided equally among the households. 

The P1 firm uses resources (RP) and labor to produce plant inventory (P1). Households, 

the H1 firm and the industrial sector are economic consumers of this P1 inventory. H2, a 

wild herbivore, preys on P1. The P1 producer devotes labor to build and maintain fences 

to keep H2 out.. The H1 firm uses P1 and P2 to produce an inventory of herbivores, H2. 

The households are the only economic consumers of this inventory. A carnivore (C1) 

preys on the H1 inventory, but the H1 firm cannot kill or otherwise interfere with C1 

because C1 is a protected species. The H1 firm’s only recourse is to invest labor in 

fences. The H1 firm pays a fee to have grazing access to P2. It is assumed that this 

access is limited, and that the H1 firm takes the maximum that it is allowed. The 

industrial sector combines resources (RP) with plants (P1) and labor to produce goods 

consumed directly by the households. Unlike the P1 and H1 inventory, which transferred 

to their respective consumers, the mass associated with the consumption of IS inventory 

is not resident in the human compartment. Rather, it goes directly to an inaccessible 

resource pool. 

The circular flow of income and spending is illustrated in Figure 2-1. Symbols used in 

this figure and throughout this report are in Appendix A. 

10
Circular Flow of Income and Spending

Firms 

P1 

H1 

IS 

Households pP1P1HH+pH1H1HH+pISISHH 

Consumption spending ($) 

Π ($) + Fixed Costs 

Entrepreneurship + Cap. Services 

Labor Market hN Labor Services d 

Wages($) = WhNd Wages($) 

Grazing fee 

Tr 

Government 

Figure 2-1. The economic model illustrating the flows of income and spending by Firms, Labor, 

Government, and Households. Note that in this model, taxes and “social costs” are not included. 

Economic Decisions of Firms and Households 

The households are the ultimate owners of the factors of production used to produce the 

goods that are traded in explicit markets. Factor incomes consist of labor incomes and 

profit incomes generated by the three industries, P1, H1 and IS. All production functions 

are specified in accordance with the law of conservation of mass: the mass emerging 

from a production process cannot exceed the total mass entering that process as inputs. 

In this initial model, the stock of physical capital is held constant. One possibility for 

introducing labor and capital into the production process is to view these inputs as simply 

reducing the amount of mass wasted in the production process. However, labor and 

capital provide a much greater service than merely reducing waste occurring in the 

production process. In particular, while labor and capital cannot create mass, they are 

nevertheless essential to the creation of value. Labor and capital provide necessary 

transformation services that alter the location of or the physical, chemical, biological, 

and/or aesthetic properties of the raw materials applied to the production function. Labor 

and capital can be substituted for each other in providing these transformation services, 

but neither labor nor capital nor a combination of the two factors can be substituted for 

raw materials in creating a given output mass. One possible exception is the extent that 

labor and capital can reduce waste, a possibility that is ignored here. 

Therefore, for this purpose, the appropriate production function views transformation 

services and raw materials as being combined in fixed-proportions. However, the 

transformation services alone may be provided by a number of combinations of labor and 

11
capital. It is further assumed that each industry and the households operate within a 

single-period planning horizon. Therefore, all saving and real investment in man-made 

capital are ignored. However, because they announce the price of their products before 

they know the actual demand for their goods, every industry may engage in unintended 

(dis)investment as they accumulate inventories (experience unfilled orders) over the 

current period. Each industry has a fixed number of equity shares (stock certificates) 

outstanding. Every household holds the same number of shares issued by a particular 

industry. Any dividend income the households receive is divided equally among the 

households. No market exists for these shares. Prices are set in terms of an abstract unit 

of account ($). Otherwise, the financial sectors are ignored. A stock of money 

presumably exists that facilitates transactions among the firms and the households, but 

any decision making as to the amount of money the agents desire to hold at any point in 

time is abstracted from. Also ignored are the social costs that the private industries 

producing plants and animals impose upon the ecosystem by growing their products, 

thereby removing nutrients from the resource pool. As a result, the market prices paid by 

the households and by the H1 industry for P1 reflect only the private out-of-pocket costs 

associated with the production of those plants. 

Also, the market price that the H1 industry charges for its product will not reflect the full 

social cost of depleting the resource pool to produce food for the domesticated 

herbivores. In addition, the households are not required to reckon with the social cost 

they impose upon the ecosystem by adding to the inaccessible resource pool as a result of 

their consumption of the various plants and animals, either through an increase in 

consumption per capita or as a result of human population growth. The H1 industry does 

pay a grazing fee to the government for access to P2, but otherwise all economic agents 

in this model ignore the ecological benefits of P2, H2 and C2. 

Households are free to decide how many hours they prefer to work based upon their 

preferences for P1, H1, IS, and leisure as well as upon the product prices, wage rate and 

the non-wage income they face. However, at the wage rate set by the IS firms, 

involuntary unemployment may result. The private industries use at least some of the 

revenue they receive from their current-period sales to pay for labor services. The P1 

industry pays part of its revenue to labor. The H1 industry also uses some of its revenue 

to pay the P1 industry for the plants it purchases from that industry and to pay a grazing 

fee to the government for access to P2. The IS industry pays for labor and for the P1 it 

buys. All three private industries then distribute any remaining profits to the households 

as dividends. In the present model, the P1, H1 and IS industries set the price of their 

respective products before trade takes place during the current period. The prices they set 

conform to their forecasts of the demand functions they will face during the current 

period. Trade takes place in the domesticated sector even though markets may fail to 

clear. Industries may hold unanticipated inventories (unfilled orders) at the end of the 

period. It is assumed that the IS industry sets the money wage rate at the beginning of the 

current period, based upon its forecasts of the household sector’s supply of labor and the 

demand for labor by the P1 and H1 industries. Involuntary unemployment (or unfilled 

job vacancies) is possible. Note that in the following equations, gi is the growth 

parameter of species i and and mi is the mortaility parameter of species i. 

12
Optimal Economic Behavior for P1 Industry 

The P1 industry applies a variable amount of labor and a fixed amount of capital to 

transform the mass of P1 that would otherwise grow naturally into a marketable product. 

The producer of P1 must also deal with the fact that H2 eats some P1 during the 

production process. The P1 industry hires labor to reduce the consumption of P1 by H2. 

The production function for P1 is given by (2-1): 

P1 = min[gP1·P1t ·RPt –mP1·P1t –cP1H2(h*

P1Nt)·H2t, P1(hP1Nt, KP1t)] (2-1) 

The production (growth) of P1 (gP1) is viewed as positively related to: (a) the size of the 

resource pool at the beginning of the period, RPt, (b) the initial stock of P1, P1t, at the 

beginning of the period, and (c) the level of transformation services, P1(·), provided by 

(variable) labor hours, hP1Nt, during the period and the (given amount of) physical 

capital, KP1

t, held by the industry at the beginning of the period. From the point of view 

of the P1 industry, the amount of P1 lost during the current period (mortality, mP1) due to 

consumption by H2, is assumed to be proportional to the initial stock of H2, H2t, with the 

size of the “consumption coefficient,” cP1H2, and negatively related to the labor hours 

devoted to reducing the amount of P1 that H2t eats during the period. Assuming the P1 

industry wastes neither the transformation services of labor and capital nor the net 

amount of P1 available for transformation (after deducting the loss due to the presence of 

H2), the amount of (transformed) P1 produced may be viewed as equal to P1(hP1Nt, KP1

t): 

P1 = P1(hP1Nt, KP1

t) (2-2) 

Note that the characteristics of P1 may differ considerably from those that would occur if 

P1 were to grow naturally without the transformation services provided by labor and 

capital. 

In addition, the amount of P1 transformed in the production process may be viewed as: 

gP1·P1t ·RPt –mP1·P1t –cP1H2(h*

P1Nt)·H2t = P1(hP1Nt, KP1t) (2-3) 

In principle (2-3) may be solved for h*

P1·Nt. To obtain the (linear approximation of) 

labor hours devoted to reducing H2’s consumption of P1 (and holding gP1, mP1 and KP1 

constant), totally differentiate (2-3) and solve for d(h*

P1·Nt ): 

d( h*

P1·Nt ) = [gP1· P1t ·d(RPt) + [gP1· RPt – mP1] ·d(P1t) –[∂P1/ ∂(hP1Nt)]·d(hP1Nt) 

– cP1H2 ·(dH2t)] / {c′P1H2 ·H2t } (2-4) 

where c′P1H2 denotes the derivative of the consumption coefficient with respect to h*

P1·Nt; 

this derivative is negative. Consequently, according to (2-5), the number of hours the P1 

industry devotes to limiting the amount of P1 consumed by H2 is negatively related to 

13
RPt and P1t, but positively related to both the number of labor hours it devotes to 

transforming P1 and the size of H2t. Therefore, the following general function is given: 

h*

P1·Nt = h*P1·Nt( P1t, RPt, H2t, hP1Nt) (2-5) 

– – + + 

The economic stock of P1 at the end of the period, P1t+1, is represented by (6): 

P1t+1 = P1t + P1(hP1Nt, KP1

t) –P1H1 –P1IS–P1HH (2-6) 

where P1H1 corresponds to the amount of P1 purchased by the H1 industry during the 

period, P1IS represents the amount of P1 purchased by the IS industry during the period 

and P1HH denotes the amount of P1 purchased by the households during the current 

period. 

Taking the current wage rate and the P1 industry’s forecast of the current period demand 

function for domesticated P1 plants as given, the P1 industry attempts to maximize the 

dividends they pay to its shareholders, subject to the restriction that the end-of-period 

stock of P1 is maintained at some predetermined level, P 1 . Then the firm’s desired 

sales of P1 (in terms of mass) is given by: 

(P1H1 +P1IS+P1HH)de = P1(hP1Nt, KP1

t) – ( P 1 –P1t ) (2-7) 

Dividends equal current net income minus net business saving. Net income is equal to 

the value of current production minus current expenses (wages and a fixed cost 

associated with the fixed stock of physical capital). Therefore, dividends are equal to the 

value of current production minus current expenses and minus net business saving. Net 

business saving is necessarily equal to the sum of the sector’s net increase in assets minus 

the net increase in liabilities during the period. In the present model, the sector’s planned 

accumulation of assets consists only of the market value of its planned accumulation of 

P1 during the current period; according to the assumptions, it does not plan to borrow 

during the current period. Therefore, current dividends are equal to the value of current 

production minus the current expenses and minus the value of the planned accumulation 

of the inventory of P1 during the period. From (2-7), the value of current production 

minus the value of the planned accumulation of inventory during the current period 

represents the value of desired current sales of P1 to other sectors during the period. 

Therefore, current dividends correspond to the current revenue from the sale of P1 to the 

producers of domesticated herbivores H1 and to the IS producers minus the current 

expenses of the P1 industry. 

Maximize: 

ΠP1 = pP1·{P1(hP1Nt, KP1

t) – ( P 1 – P1t) } 

 –W·[ hP1·Nt + h*

P1·Nt(P1t, RPt, H2t, hP1Nt)] –FP1 (2-8) 

where FP1 denotes the current period’s fixed cost. Since the P1 industry is viewed as 

14
owning the P1 that it produces and harvests, the amount of P1 harvested does not enter 

directly into its cost considerations because the industry does not voluntarily consider the 

social cost of its production of P1. 

At the beginning of the period, the P1 industry announces its product price, pP1. Let the 

P1 industry’s forecast of the current period demand for its product be represented by the 

following demand function: 

P1de = P1de( pP1, ωH1, ωIS, ωP1HH) = P1H1de(pP1, ωH1) + P1ISde(pP1, ωIS) (2-9) 

+(P1HH/N)de(pP1, ωP1HH)·Nt 

This demand function denotes the P1 industry’s forecasted total (market) demand for P1 

by the H1 and IS industries and by the households. It is assumed that the derivative of 

P1de with respect to pP1 is negative (the sign will be verified later when the economic 

behavior of the H1 and IS industries and the households is specified); ωH1 , ωIS and ωP1HH 

denote vectors of shift parameters. Solving the inverse function for pP1 yields: 

pP1 = pP1(P1de, ωH1, ωIS, ωP1HH, Nt ) (2-10) 

– + + + + 

where: 

∂pP1/∂(P1de) ≡ 1/ [∂(P1de)/ ∂pP1 ] < 0 

∂pP1/∂ωi ≡ – [∂(P1de)/ ∂ωi ]/[∂(P1de)/ ∂pP1 ] > 0 

and ∂pP1/∂Nt ≡ – [(P1HH/N)de]/[∂(P1de)/ ∂pP1 ] > 0. (2-11) 

An increase in ωi or Nt represents an outward shift in the P1 industry’s forecast of the 

market demand for its product. 

The price pP1 given by (2-10) represents the maximum uniform price the P1 sector 

expects it can value each alternative quantity of P1 it produces during the period. Based 

upon this function, the sector’s forecasted current revenue function is given by: 

RP1e = pP1· P1de = pP1(P1de, ωH1, ωIS, ωP1HH, Nt ) · P1de

 = RP1e(P1de, ωH1, ωIS, ωP1HH, Nt ) (2-12) 

Substituting the right hand side of (2-7) for P1de in (2-12) yields: 

RP1e

 = RP1e[P1(hP1Nt, KP1

t) –( P 1– P1t), ωH1, ωIS, ωP1HH, Nt] (2-13) 

Substituting the right hand side of (2-13) for the first term on the right hand side of (2-8) 

yields: 

ΠP1 = RP1e

{P1(hP1Nt, KP1

t) –( P 1 – P1t), ωH1, ωIS, ωP1HH, Nt } 

 –W·[hP1·Nt + h*

P1·Nt(P1t, RPt, H2t, hP1Nt)] –FP1 (2-14) 

15
 

The only independent choice variable for the industry in (2-14) is the number of hours 

that its employees work transforming P1, hP1Nt. Maximizing (2-14) with respect to 

hP1·Nt yields the necessary condition: 

[∂RP1e

/∂(P1de)] ·[ ∂(P1)/∂(hP1Nt)] = W·[ 1 + [ ∂(h*

P1·Nt)/∂(P1)]· [∂(P1)/∂(hP1·Nt)]] 

= W·[ 1 – [∂(P1)/∂(hP1·Nt)]/{c′P1H2 ·H2t }] (2-15) 

The left hand side of (2-15) represents the marginal revenue product of labor used in the 

P1 industry to transform P1 into a marketable product, and the condition requires the 

profit-maximizing firms in the industry to hire labor up to the point at which the marginal 

revenue product of labor equals the money wage times Equation 2-1 plus the extra labor 

required to reduce the amount that H2 eats as the sector adds to P1. Holding P1t, RPt, 

and H2t constant, if the firm hires more labor to transform P1 into a marketable product, 

it must also hire more labor to reduce the loss of P1 due to its consumption by H2, in 

order to obtain the extra P1 that will be used as the raw material input. 

Assuming a diminishing marginal product of labor in transforming P1, a diminishing 

marginal product of labor in reducing loss, and that marginal revenue decreases as the 

number of units sold increases, then from (2-15), the profit maximizing level of labor to 

be used in the P1 industry becomes a decreasing function of the money wage. The P1 

industry’s demand functions for both types of labor are given by (2-16) and (2-17): 

(hP1·Nt)

d

 = hP1·Nt

d

(W, ( P 1 – P1t), ωH1, ωIS, ωP1HH, Nt, P1t, RPt, H2t) (2-16) 

– + + + + + + + – 

(h*

P1·Nt)

d

 = h*

P1·Nt

d

(W, ( P 1 – P1t), ωH1, ωIS, ωP1HH, Nt, P1t, RPt, H2t) (2-17) 

– + + + + + – – + 

From (2-16) an increase in ( P 1 – P1t) reduces, ceteris paribus, the amount of P1 the 

sector plans to sell, thereby raising the marginal revenue product of labor used to 

transform P1; an increase in ωH1, ωIS, ωP1HH, or Nt shifts the market demand for P1 

outward, thereby increasing the price that buyers are willing to pay for a given amount of 

P1, and raising the forecasted marginal revenue product for a given amount of labor used 

to transform P1. The effect upon the amount of labor used to conserve P1 responds in the 

same direction as the demand for labor used to transform P1. However, holding the 

amount of labor used to transform P1 unchanged, in accordance with (2-3), an increase in 

either P1t or RPt reduces, ceteris paribus, the amount of labor required to conserve P1, 

causing the demand for that type of labor to diminish. As the quantity of h*

P1·Nt 

decreases, its marginal product increases, thereby reducing the marginal cost of adding a 

unit of hP1·Nt; the industry’s demand for hP1·Nt increases with an increase in P1t or RPt. 

The larger the stock of H2t at the beginning of the period, however, the greater will be the 

industry’s demand for h*

P1·Nt and the smaller will be its demand for hP1·Nt. Substituting 

the right hand side of (2-16) into the P1 industry’s transformation function, (2-2), yields 

the P1 industry’s current production of domesticated plants consistent with profit 

maximization in that industry: 

16
P1s = P1s (hP1·Nt

d

[W, ( P 1 – P1t), ωH1, ωIS, ωP1HH, Nt, P1t, RPt, H2t], KP1

t) (2-18) 

The price that the P1 industry announces for the current period is then given by: 

pP1 = pP1{P1de, ωH1, ωIS, ωP1HH, Nt}

 = pP1{P1s (hP1·Nt

d

[W, ( P 1 – P1t), ωH1, ωIS, ωP1HH, Nt, P1t, RPt, H2t], KP1

t) 

–( P 1– P1t), ωH1, ωIS, ωP1HH, Nt} (2-19) 

The actual inventory of P1 at the end of the period is given by: 

P1t+1 = P1t + P1s –P1H1 –P1IS –P1HH. (2-20) 

Optimal Economic Behavior for H1 Industry 

H1 industry buys P1, labor and grazing rights to P2 and sells output to the households. 

The assumptions are that the government sector fixes the grazing fee, p2, and that the H1 

industry buys the maximum amount of P2, denoted by P2H1, permitted by the 

government. C1 also consumes H1; the H1 industry hires labor to limit the amount of H1 

eaten by C1. 

The production function for H1 is shown by (2-21): 

H1 = min[P1H1 + P2H1 –mH1·H1t –cH1C1(h*

H1Nt) ·C1t , H1(hH1Nt, KH1t)] (2-21) 

where H1(·) denotes the transformation services provided by labor and capital. 

Assuming that the H1 industry wastes neither the transformation services provided by 

labor and capital nor the net amount of H1 available for transformation, the amount of H1 

available to the market during the current period may be written as: 

H1 = H1(hH1Nt, KH1

t) (2-22) 

and also: 

 P1H1 + P2H1 –mH1·H1t –cH1C1(h*

H1Nt) ·C1t = H1(hH1Nt, KH1t) (2-23) 

The consumption coefficient, cH1C1, is positive; the first derivative of this coefficient is 

negative and its second derivative is positive. 

In principle, (2-23) may be solved for h*

H1Nt in terms of hH1Nt, P1H1, C1t , and P2H1. 

Totally differentiating (2-23) yields the following: 

d(h*

H1Nt) = [d(P1H1) +d(P2H1) –mH1·d(H1t) –[∂H1/ ∂(hH1Nt)]·d(hH1Nt) 

– cH1C1 ·d(C1t)] / {c′H1C1 ·C1t } (2-24) 

17
Therefore, in general: 

h*

H1Nt = h*H1Nt(P1H1, P2H1, H1t, hH1Nt, C1t) (2-25) 

– – + + + 

The stock of H1 at the end of the current period is then given by: 

H1t+1 = H1t + H1(hH1Nt, KH1

t) –H1HH. (2-26) 

H1HH denotes the amount of H1 sold to the households during the period. 

Taking as given the current wage rate and its forecast of the current period demand 

function for domesticated herbivores, H1, the H1 firms attempt to maximize the 

dividends to their shareholders, subject to the restriction that they maintain the end-ofperiod stock of H1 at some predetermined level, H 1 . The firm’s desired sales of H1 (in 

terms of mass) is given by: 

(H1HH)de = H1(hH1Nt, KP2

t) –( H 1 –H1t ) (2-27) 

Therefore, the objective of the H1 industry is to maximize: 

ΠH1 = pH1·{H1(hH1Nt, KH1

t) – ( H 1 – H1t) } –pP1·P1H1 –pP2·P2H1

 –W·[hH1Nt + h*

H1Nt(P1H1, P2H1, H1t, hH1Nt, C1t)] –FH1 (2-28) 

where FH1 represents the industry’s fixed cost. 

The product price, pH1, is announced by the H1 industry at the beginning of the period. 

Let the H1 industry’s forecast of the current period demand for its product be given by 

the following demand function: 

H1HHde = (H1HH/Nt)

de(pH1, ωH1HH)·Nt (2-29) 

This demand function denotes the H1 industry’s forecast of the households’ demand for 

H1. The assumption is that the derivative of H1de with respect to pH1 is negative (the sign 

will be verified later when the economic behavior of the households is specified); ωH1HH 

denotes a vector of shift parameters. Solving the inverse function for pH1 yields: 

pH1= pH1(H1HHde, ωH1HH, Nt) (2-30) 

where 

∂pH1/∂(H1HHde) ≡ 1/ {[∂(H1HH/N)de/∂pH1]·N} <0 (2-31) 

∂pH1/∂ωH1HH ≡ – [∂(H1HH/N)de/ ∂ωH1HH ]/{[∂(H1HHde)/ ∂pH1 ]·N} > 0 (2-32) 

∂pH1/∂Nt ≡ – [(H1HH/N)de]/{[∂(H1HHde)/ ∂pH1 ]·N} > 0 (2-33) 

and where a rise in ωH1HH or Nt represents an outward shift in the demand for H1 function 

18
facing the H1 industry. 

The price pH1 given by (2-30) represents the maximum uniform price the H1 sector 

expects it can value each alternative quantity of H1 it produces during the period. Based 

upon this function, the sector’s forecasted current revenue function is given by: 

e RH1 = pH1·H1HHde = pH1(H1HHde, ωH1HH, Nt)·H1HHde = RH1e

(H1HHde, ωH1HH, Nt ) (2-34) 

Substituting the right hand side of (2-27) for H1HHde in (2-34) yields: 

RH1e

 = RH1e[H1(hH1Nt, KP2

t) –( H 1 –H1t ), ωH1HH, Nt ] (2-35) 

Substituting the right hand side of (2-35) for the first term on the right-hand-side of (2­

28) yields: 

ΠH1 = RH1e{H1(hH1Nt, KP2

t) –( H 1 –H1t ), ωH1HH, Nt} –pP1·P1H1 –pP2·P2H1 

 –W·[hH1Nt + h*

H1Nt(P1H1, P2H1, H1t, hH1Nt, C1t)] –FH1 (2-36) 

Given H 1 , two choice variables confront the H1 industry in this simplified model: the 

number of hours that employees work in the industry transforming P1 and P2 into H1; 

and the amount of the intermediate good, P1, the H1 industry buys from the P1 industry. 

The amount of labor the H1 industry employs to reduce C1’s consumption of H1 then 

follows from (2-25). Maximizing (2-36) with respect to hH1·Nt yields the necessary 

condition: 

[∂RH1e/∂(H1HHde)]·[∂(H1)/∂(hH1Nt)] =W·[1+[∂(h*

H1Nt)/∂(H1)]·[∂(H1)/∂(hH1Nt)]]

 = W·[1–{[∂(H1)/∂(hH1Nt)]/(c′H1C1 ·C1t )}] (2-37) 

The left hand side of (2-37) represents the marginal revenue product of labor used in the 

H1 industry to transform mass into marketable H1. The condition requires the profitmaximizing firms in the industry to hire this type of labor up to the point at which the 

marginal revenue product of labor equals the money wage, plus the marginal cost of 

adding as well the requisite amount of labor to further limit the consumption of H1 by C1 

in order to obtain, ceteris paribus, the extra H1 mass to be transformed. 

The first-order condition for P1H1 is given by: 

–W·[ ∂ (h*

H1Nt)/ ∂(P1H1)] = pP1 (2-38) 

or: 

 –W/{c′H1C1 ·C1t } = pP1 (2-39) 

According to condition (2-38), the industry should continue to buy P1 up to the point at 

which the marginal revenue to the industry from an extra unit of P1, namely the wages it 

saves because it can reduce its use of labor to limit C1’s consumption of H1, is equal to 

the price of p1. 

19
Assuming a diminishing marginal product of labor and assuming that marginal revenue 

decreases as the number of units sold increases, then taking total differentials of (2-37) 

and (2-39) yields a set of equations that can be solved for the responses in the choice 

variables to various parametric changes. In particular, the following demand functions 

are found for labor hours, hH1Nt, and for P1H1: 

(hH1·Nt)

d

 = hH1·Nt

d

(W, pP1 , ( H 1 –H1t ), ωH1HH, Nt) (2-40) 

– – + + 

(P1H1)d = P1H1d(W, pP1, ( H 1 –H1t ), ωH1HH, Nt) (2-41) 

– – + + 

Since, ceteris paribus, the inputs P1H1 and labor that is used to transform H1, hH1Nt, are 

technical complements, an increase in the price of either input results in the H1 industry 

demanding less of both inputs. Substituting the right-hand sides of (2-40) and (2-41) into 

(2-25) yields the H1 industry’s demand for labor for the purpose of limiting the amount 

of H1 consumed by C1: 

(h*

H1Nt)

d = (h*

H1Nt)

d [(P1H1)d

(W, pP1, ( H 1 –H1t ), ωH1HH, Nt ), P2H1, H1t, 

(2-42) 

(hH1Nt)

d

(W, pP1, ( H 1 –H1t ), ωH1HH, Nt) , C1t] 

or: (2-43) 

(h*

H1Nt)

d = (h*

H1Nt)

d [W, pP1, ( H 1 –H1t ), ωH1HH, Nt P2H1, H1t, C1t] 

– + ± ± ± – + + 

where it has been assumed that the own-price effects upon (P1H1)d and (hH1Nt)

d 

dominate the indirect, or cross-price effects. 

Substituting the right hand side of (2-40) into (2-22) yields the amount of H1 the industry 

desires to supply to the market during the current period: 

H1s = H1s

[hH1·Nt

d

(W, pP1 , ( H 1 –H1t ), ωH1HH, Nt) , KH1

t] 

– – + + + 

= H1s(W, pP1 , ( H 1 –H1t ), ωH1HH, Nt) , KH1

t) (2-44) 

– – + + 

The price that the H1 industry announces for the current period is then given by: 

pH1= pH1(H1HHde, ωH1HH) 

= pH1(H1s

(W, pP1 , ( H 1 –H1t ), ωH1HH, Nt) , KH1

t) –( H 1 –H1t ), ωH1HH, Nt) (2-45) 

– – + + 

The actual inventory of H1 at the end of the period is given by: 

H1t+1 = H1t + H1s –H1HH. (2-46) 

20
Optimal Economic Behavior for IS Industry 

The IS industry buys P1 and combines P1 and RP in fixed proportions. Since RP is 

“free” the industry would use only RP if it could produce IS using variable proportions of 

P1 and RP. P1 and RP are combined with a variable amount of labor that is used to 

transform P1 and RP into IS. This arrangement is slightly more complicated than the H1 

industry, since both P1 and RP vary. The IS industry buys P1, combines it with RP using 

variable labor and a fixed amount of capital to produce IS, which it sells to the 

households. Therefore the production function for IS may be written as: 

IS = min[P1IS/θ, RP/λ, IS(hISNt, KIS

t)] (2-47) 

where θ denotes the amount of P1 necessary to produce a unit of IS and λ denotes the 

amount of RP necessary to produce a unit of IS. The function IS(·) denotes the 

transformation services provided by labor and capital in the process of producing a 

marketable unit of IS. 

Assuming that the IS industry wastes neither the transformation services provided by 

labor and capital, nor P1 nor RP, the number of units of IS available to the market during 

the current period may be written as: 

IS = IS(hISNt, KIS

t) (2-48) 

where, 

P1IS = θ·IS(hISNt, KIS

t) (2-49) 

RPIS = λ·IS(hISNt, KIS

t) (2-50) 

In terms of mass, the stock of IS at the end of the current period is given by: 

(θ+λ)ISt+1 = (θ+λ)ISt +(θ+λ)·IS(hISNt, KIS

t) –(θ+λ)·ISHH. (2-51) 

where ISHH denotes the number of units of IS purchased by the households and 

(θ+λ)·ISHH represents the amount of IS purchased by the households during the period in 

terms of mass. 

Assuming that the IS industry sets the current wage rate, the assumption is that the 

industry formulates an anticipated net supply of labor function, NS

, which presumably 

reflects the industry’s estimate of (a) the amount of labor the household sector is willing 

to supply to the labor market during the current period at each money wage, W, less (b) 

the amount of labor that the P1 and H1 industries demand at each money wage. It is 

assumed that the estimated net supply function is given by: 

21
(hISN)se = (hISN)se(W) (2-52) 

+ 

Taking the industry’s estimated labor supply function and its forecast of the household 

sector’s current period demand function for IS as given, the IS firms attempt to maximize 

the dividends to their shareholders, subject to the restriction that they maintain the endof-period stock of IS at some predetermined level, I S . The firm’s desired sales of IS 

(in terms of units of IS) is given by: 

(ISHH)de

 = IS[(hISN)se(W) , KIS

t] –( I S –ISt ) (2-53) 

Therefore, the objective of the IS industry is to maximize: 

ΠIS = pIS·{IS[(hISN)se(W) , KIS

t] –( I S –ISt )} –W· (hISN)se(W) –FIS –pP1·P1IS (2-54) 

or: 

ΠIS = pIS·{IS[(hISN)se(W) , KIS

t] –( I S –ISt )} 

 –W· (hISN)se(W) –FIS –pP1·θ· IS[(hISN)se(W) , KIS

t] (2-55) 

where FIS represents the industry’s fixed cost. 

The product price, pIS, is announced by the IS industry at the beginning of the period. Let 

the IS industry’s forecast of the current period demand for its product be given by the 

following demand function: 

ISHHde = (ISHH/N)de(pIS, ωISHH)·Nt (2-56) 

– + 

This demand function denotes the IS industry’s forecast of the households’ demand for 

IS. It is assumed (it will be verified later when the economic behavior of the households 

is specified) that the derivative of (ISHH/N)de with respect to pIS is negative; ωISHH 

denotes a vector of shift parameters. Solving the inverse function for pIS yields: 

pIS= pIS(ISHHde, ωISHH, Nt ) (2-57) 

– + + 

where: 

∂pIS/∂(ISHHde) ≡ 1/{[∂(ISHH/N)de)/∂pIS ]·Nt} < 0 

∂pIS/∂ωISHH ≡ – [∂(ISHH/N)de/∂ωISHH ]/{[∂(ISHH/N)de/∂pIS ]·Nt} > 0 

and 

∂pIS/∂Nt ≡ – (ISHH/N)de/{[∂(ISHHde)/∂pIS ]·Nt} > 0 

and where a rise in ωISHH or Nt represents an outward shift in the demand for IS function 

22
facing the IS industry. 

The price pIS given by (2-57) represents the maximum uniform price the IS sector expects 

it can value each alternative amount of IS (in units) it produces during the period. Based 

upon this function, the sector’s forecasted current revenue function is given by: 

e RIS = pIS· ISHHde = pIS(ISHHde, ωISHH, Nt)·ISHHde = RISe

(ISHHde, ωISHH, Nt) (2-58) 

Substituting the right hand side of (53) for ISHHde in (58) yields: 

RISe

 = RISe[IS[(hISN)se(W) , KIS

t] –( I S –ISt ), ωISHH, Nt] (2-59) 

Substituting the right-hand-side of (2-59) for the first term on the right hand side of 

(2-55) yields: 

ΠIS = RISe[IS[(hISN)se(W) , KIS

t] –( I S –ISt ), ωISHH, Nt ] 

–W· (hISN)se(W) –FIS –pP1·θ· IS[(hISN)se(W) , KIS

t] (2-60) 

Given I S , the only choice variable for the IS industry in this simplified model consists 

of the wage rate that it announces to attract the appropriate number of employees to 

transform P1 and RP into IS products; the amount of the intermediate good P1 the IS 

sector buys, P1IS, must then conform to the amount of mass necessary to produce that 

amount of IS. Maximizing (2-60) with respect to W yields the necessary condition: 

[∂RISe/∂(ISHHde)]·[∂(IS)/∂(hISNt)]·[∂(hISNt)/∂W]

 = (hISN)se(W) +W·[∂(hISNt)/∂W]+ pP1·θ·[∂(IS)/∂(hISNt)]·[∂(hISNt)/∂W] (2-61) 

The left hand side of (2-61) represents the marginal revenue product of labor in the IS 

industry in terms of transforming P1 and RP into IS; the condition requires that the profitmaximizing firms in the industry continue to raise the wage rate, and therefore continue 

to hire labor up to the point at which the marginal revenue product of from the extra labor 

induced by the higher wage rate equals the marginal cost of raising the money wage, 

(hISN)se(W) +W·[∂(hISNt)/∂W], plus the marginal cost of adding as well the requisite 

amount of P1 to the production process to obtain the extra IS. The extra amount of RP 

that must accompany the extra P1 is considered to be a free good and is given by (2-50). 

Assuming a diminishing marginal product of labor and that marginal revenue decreases 

as the number of units sold increases, then from (2-61), the profit maximizing wage to be 

announced by the IS industry is a decreasing function of the price of pP1 and an 

increasing function of the shift parameters, ( I S –ISt ) and ωISHH. The money wage 

function for the IS industry becomes: 

 W = W( pP1·θ, ( I S –ISt ), ωISHH, Nt ) (2-62) 

– + + + 

Substituting the right hand side of (2-62) into (2-53) yields the amount of IS the industry 

23
desires to supply to the market during the current period: 

(ISHH)s = IS[(hISN)se[W (pP1·θ, ( I S –ISt ), ωISHH)] , KIS

t] –( I S –ISt ) (2-63) 

Substituting the right hand side of (2-63) into (2-57) yields the price pIS announced by the 

IS industry: 

pIS= pIS(IS[(hISN)se[W (pP1·θ, ( I S –ISt ), ωISHH, Nt )] , KIS

t], ωISHH, Nt ) (2-64) 

In terms of mass, the stock of IS at the end of the period is given by: 

(θ+λ)·ISt+1 = (θ+λ)·ISt + (θ+λ)·IS[(hISN)se[W(pP1·θ, ( I S –ISt ), ωISHH, Nt )] , KIS

t] 

–(θ+λ)·ISHH. (2-65) 

From (2-49) and (2-50), the amount of mass transferred from P1 and RP to IS is then 

given by: 

P1IS = θ·IS((hISNt)

d

 (W, KIS

t, pP1·θ, ( I S –ISt ), ωISHH, Nt ), KISt) (2-66) 

RPIS = λ·IS((hISNt)

d

 (W, KIS

t, pP1·θ, ( I S –ISt ), ωISHH, Nt ), KISt) (2-67) 

Government Sector 

The only revenue the government sector receives is from the grazing rights it grants to 

the H1 industry equal to pP2·P2H1. For the time being, the assumption is that the 

government transfers this revenue to the households: 

Tr = pP2·P2H1. (2-68) 

GDP 

The level of nominal GDP is defined as the total market value of all units of final goods 

produced in the economy during the period, which also corresponds to the sum of all 

spending on final goods plus the values of the changes in inventories of all goods: 

GDP = pP1·{P1HH} +pH1·{H1HH} + pIS·{ISHH} 

+ pP1(P1t+1 –P1t) +pH1(H1t+1 –Ht) +pIS(ISt+1 –ISt) (2-69) 

Optimal Economic Behavior of the Household Sector 

In the present section the household sector’s market demands for P1, H1, and IS as well 

as its market supply of labor is obtained. It is assumed that the objective of the household 

sector is to maximize current utility, which is assumed to be a positive increasing 

24
function of consumption per capita of P1, H1, IS and leisure, l. It is assumed that 

marginal utility is positive, but decreasing. For simplicity, cross-partials are ignored. 

Maximize: 

U = U(P1HH/Nt , H1HH/Nt , ISHH/Nt, l ) (2-70) 

The households attempt to maximize utility subject to their budget constraint. This 

constraint stipulates that total household spending on all three products equals its current 

income, consisting of wage and non-wage income. It is assumed that the households take 

their current non-wage incomes as given, but attempt to choose their wage income. The 

households are not permitted to make adjustments based on how their current non-wage 

income may be affected by their work/leisure choice, or by the market basket of goods 

they decide to purchase. In terms of the abstract unit of account, the budget constraint 

facing the household sector is given by: 

Π +Tr +W·(h*–l)·Nt = pP1·(P1HH/Nt )·Nt + pH1·(H1HH/Nt )·Nt +pIS·(ISHH/Nt)·Nt (2-71) 

where: 

Π = Π PP1 +FPP1 + Π H1 +FH1 +Π IS +FIS. (2-72) 

Equation (2-70) is rewritten by dividing all terms by WNt, thereby expressing the terms 

in the constraint in units of labor per capita: 

(Π +Tr)/(W·Nt) +(h*–l) = (pP1/W)·(P1HH/Nt )+ (pH1 /W)·(H1HH/Nt ) +(pIS/W)·(ISHH/Nt). (2-73) 

The following Lagrangean function is set up: 

Max U = U(P1HH/Nt , H1HH/Nt , ISHH/Nt, l ) 

+λ{(Π+Tr) /(W·Nt) +(h*

–l) –(pP1/W)·(P1HH/Nt ) –(pH1 /W)·(H1HH/Nt ) 

–(pIS/W)·(ISHH/Nt)} (2-74) 

The first-order necessary conditions for the maximization of (2-74) with respect to 

P1HH/Nt, H1HH/Nt, ISHH/Nt and l are given by (2-75)–(2-78). From (2-71), the sum 

inside braces in (2-74) must equal zero; this restriction is added as condition (2-79): 

∂U/∂(P1HH/Nt) –λ (pP1/W) = 0. (2-75) 

∂U/∂(H1HH/Nt) –λ (pH1/W) = 0. (2-76) 

∂U/∂(ISHH/Nt) –λ (pIS/W) = 0. (2-77) 

∂U/∂l –λ = 0. (2-78) 

(Π+Tr)/(W·Nt) +(h*–l) –(pP1/W)·(P1HH/Nt ) 

–(pH1 /W)·(H1HH/Nt ) –(pIS/W)·(ISHH/Nt ) = 0. (2-79) 

According to (2-75)–(2-79), only non-wage income (in labor units, (Π+Tr)/W), and the 

relative prices, pP1/W, pH1/W, and pIS/W (expressed here in labor units as well) matter in 

25
terms of the optimal combinations of l, P1HH/Nt, H1HH/Nt and ISHH/Nt selected by the 

households. The budget constraint, expression (2-79), is homogeneous of degree zero in 

W, Π+Tr, pP1, pH1 and pIS. A doubling of all prices (including W) and non-wage income 

is consistent with the same combination of l, P1HH/Nt , H1HH/Nt and ISHH/Nt. Taking 

the total differentials of each of (2-75)–(2-79) yields the following system of equations, 

which are written in matrix form as (2-80): 

=

d(P1HH/Nt u ) 11 u12 u13 u14 –(pP1/W)

d(H1HH/Nt) u21 u22 u23 u24 –(pH1/W) 

d(ISHH/Nt) u31 u32 u33 u34 –(pIS/W) 

d(l ) u41 u42 u43 u44 –1 

dλ –(pP1/W) –(pH1/W) –(pIS/W) –1 0 

λ d(pP1/W ) 

λ d( pH1/W) 

λ d( pIS/W) 

0 

–d[(Π+Tr)/(W·Nt)] +(P1HH/Nt )·d(pP1/W) +(H1HH/Nt )·d(pH1/W) +(ISHH/Nt )·d(pIS/W) 

(2-80) 

The coefficients uij denote the second partials of utility; the own-second partials, uii, are 

assumed to be negative. For simplicity, the cross-partials are treated as small relative to 

the own-partials as the signs of the effects of the parametric changes upon the choice 

variables are determined. Pre-multiplying the column vector on the right hand side of (2­

79) by the inverse of the bordered Hessian matrix yields total differentials of the 

household sector’s demands for P1HH/Nt, H1HH/Nt, ISHH/Nt and l. Assuming that the 

substitution effects dominate the income effects, these total differentials indicate the 

following responses in the household sector’s demand functions for P1HH/Nt, H1HH/Nt, 

ISHH/Nt and l respectively: 

(P1HH/Nt)

d

 = (P1HH/N)d[(Π+Tr)/WNt, pP1/W, pH1/W, pIS/W ) (2-81) 

+ – + + 

(H1HH/Nt)

d

 = (H1HH/N)d(Π+Tr/WNt, pP1/ W, pH1/W, pIS/W) (2-82) 

+ + – + 

26
(ISHH/Nt)

d

 = (ISHH/N)d(Π+Tr/WNt, pP1/ W, pH1/W, pIS/W) (2-83) 

+ + + – 

l

d

 = l

d

( Π/WNt, pP1/W, pH1/W, pIS/W) (2-84) 

+ + + + 

Expression (2-79) imposes restrictions on the responses in (2-81)–(2-84). From (2-79), 

and the fact that the hours spent working, h, equals total hours in the period, h*, minus 

leisure time, l, the household sector’s per capita supply of labor, hs, is given by (2-85): 

hs = hs[ (Π+Tr)/WNt, pP1/W, pH1/W, pIS/W). (2-85) 

– – – – 

Substituting (2-81)–(2-84) into constraint (2-79) yields: 

(Π+Tr)/(W·Nt) +h* – l

d

((Π+Tr)/WNt, pP1/W, pH1/W, pIS/W)

 – (pP1/W)·(P1HH/N)d((Π+Tr)/WNt, pP1/W, pH1/W, pIS/W)

 –(pH1/W)·(H1HH/N)d((Π+Tr)/WNt, pP1/ W, pH1/W, pIS/W)

 –(pIS/W)·(ISHH/N)d((Π+Tr)/WNt, pP1/ W, pH1/W, pIS/W) = 0. (2-86) 

Since (2-85) must hold identically, the following restrictions apply: 

1 ≡ ∂l

d

/∂(Π/WNt) +(pP1/W)·[∂(P1HH/N)d/∂(Π/WNt)] 

+(pH1/W)·[∂(H1HH/N)d/∂(Π/WNt)] +(pIS/W)·[∂(ISHH/N)d/∂(Π/WNt)] (2-87) 

0 ≡ ∂l

d

/∂(pP1/W) +(pP1/W)·[∂(P1HH/N)d/∂(pP1/W)] +(P1HH)

d

 +(pH1/W)·[∂(H1HH/N)d/∂(pP1/W) +(pIS/W)·[∂(ISHH/N)d/∂(pP1/W)] (2-88) 

0 ≡ ∂l

d

/∂(pH1/W) +(pP1/W)·[∂(P1HH/N)d/∂(pH1/W)]

 +(pH1/W)·[∂(H1HH/N)d/∂(pH1/W) +(H1HH)d +(pIS/W)·[∂(ISHH/N)d/∂(pH1/W)] (2-89) 

0 ≡ ∂l

d

/∂(pIS/W) +(pP1/W)·[∂(P1HH/N)d/∂(pIS/W)]

 +(pH1/W)·[∂(H1HH/N)d/∂(pIS/W) +(pIS/W)·[∂(ISHH/N)d/∂(pIS/W)] +(ISHH)d (2-90) 

Gross Domestic Product 

As shown in expression (2-69), repeated here for convenience, nominal gross domestic 

product, GDP, is defined as the total market value of all units of final goods and services 

produced during a period of time. In the present context, the market value of the goods 

the households purchase from the private firms plus the changes in inventories of all 

goods represent nominal GDP: 

GDP = pP1·{P1HH} +pH1·{H1HH} + pIS·{ISHH} 

+ pP1(P1t+1 –P1t) +pH1(H1t+1 –Ht) +pIS(ISt+1 –ISt) (2-69) 

Profits are defined as value added in production minus current expenses. In this model, 

current expenses consist of wages, rental income to the owners of physical capital (“fixed 

cost”), and fees paid to the government for grazing rights. Therefore, profits plus wages, 

plus rental income to the owners of physical capital, and fees paid to the government for 

27
grazing rights are equal to the value added in production (value of the goods produced in 

a stage of production minus the purchase of intermediate goods entering that stage). In 

this model, the fees paid by the H1 industry for grazing rights are returned to the 

households as transfer payments. Therefore, the sum of profits plus wages plus rental 

incomes plus transfers is equal to the value of production minus the purchase of 

intermediate goods. Since the value of production coincides with the value of sales plus 

the change in inventories, the following three relationships hold for the P1, H1, and IS 

industries respectively: 

Π*

P1 +FP1 +W·[hP1·Nt +h*P1·Nt] = pP1·{P1HH +P1H1 +P1IS + (P1t+1 –P1t)} (2-91) 

Π*

H1 +FH1 +W·[hH1Nt +h*H1Nt] +Tr = pH1·{H1HH +(H1t+1 –Ht)} –pP1·P1H1 (2-92) 

Π*

IS +FIS +W·(hISN) = pIS·{ISHH +( ISt+1 –ISt )} –pP1·P1IS (2-93) 

Summing across (91)–(93) yields: 

Π** + Tr +WhNt = pP1·{P1HH + (P1t+1 –P1t) }+ pH1·{H1HH +(H1t+1 –Ht) } 

+pIS·{ISHH +(ISt+1 –ISt )} (2-94) 

Since the right hand side of (2-94) corresponds to the right hand side of (2-69), nominal 

GDP is also equal to total non-wage income, Π**, plus transfers plus total wages: 

GDP = Π** +Tr +W·h·N. (2-95) 

Note that non-wage income corresponds to profits plus the rental income to owners: 

Π** = Π*

P1 +Π*H1 +Π*IS +FP1 +FH1 +FIS (2-96) 

Profits of the individual industries correspond to the dividends they pay to the households 

plus net business saving (additions to retained earnings). In this model net business 

saving corresponds to the change in inventories. Therefore, 

Π*

P1 = ΠP1 + pP1·(P1t+1 –P1t) (2-97) 

Π*

H1 = ΠH1 + pH1·(H1t+1 –H1t) (2-98) 

Π*

IS = ΠIS + pIS·(ISt+1 –ISt) (2-99) 

Real GDP, or RGDP, written in terms of units of goods, is given by: 

RGDP = P1HH + H1HH +ISHH +(P1t+1 –P1t) +(H1t+1 –Ht) +(ISt+1 –ISt ) (2-100) 

Naturally Occurring Changes in Biological and Physical Resources 

Growth of Non-domesticated Plants, P2 

The production (natural growth) of non-domesticated plant P2 during the current period 

28
is viewed as proportional to both the volume of P2 at the beginning of the current period 

and the size of the resource pool at the beginning of the current period. In addition, since 

plant P2 is able to draw mass from the inaccessible resource pool, plant P2’s growth is 

also proportional to the size of the inaccessible resource pool as well. But the growth of 

P2 is diminished by the amount of the plant eaten by the domesticated, H1, and nondomesticated herbivores, H2 and H3, during the period. Net growth of P2 is shown by 

(2-101). 

P2 = P2t·(gP2·RPt +gIRPP2·IRPt) –mP2·P2t –P2H1 –H2t·gH2·P2t –H3t·gP2H2·P2t. (2-101) 

The volume of P2 at the end of the current period is equal to the volume at the beginning 

of the current period, plus the net natural growth of P2, minus the amount of the plant 

eaten by the domesticated and non-domesticated herbivores, H2 and H3, during the 

period: 

P2t+1 = P2t +P2t·(gP2·RPt +gIRPP2·IRPt) –mP2·P2t –P2H1 –H2t·gH2·P2t –H3t·gP2H2·P2t. (2-102) 

Growth of Non-domesticated Herbivores, H2 

The production (natural growth) of H2 during the current period is viewed as 

proportional to both the amount of H2 at the beginning of the current period and the 

amount of plant P2 at the beginning of the current period. In addition the growth of H2 is 

increased by the amount of P1 that H2 eats during the current period, but diminished by 

the amount of H2 eaten by C1 and C2. Net growth is given by (2-103): 

H2 = P1H2 +H2t·gH2·P2t –mH2·H2t –H2C1 –C2t·gH3C2·H3t (2-103) 

and the end-of-period stock of H2 is given by: 

H2t+1 = H2t + P1H2 +H2t·gH2·P2t –mH2·H2t –H2C1 –C2t·gH3C2·H3t. (2-104) 

Growth of Non-domesticated Carnivores, C1 

The mass of C1 increases during the period by the amount of mass that C1 eats of H1, 

cH1C1·C1t, and H2, H2C1; its mass decreases by its natural mortality, mC1·C1t. Net 

growth of C1 during the period is represented by (2-105) and the end-of-period stock of 

C1 is given by (2-106): 

C1 = cH1C1·C1t + H2C1 –mC1·C1t (2-105) 

and the end-of-period stock of C1 is given by: 

C1t+1 = C1t + cH1C1·C1t + H2C1 –mC1·C1t (2-106) 

29
Growth of Non-domesticated Plants, P3 

The production (natural growth) of non-domesticated plant P3 during the current period 

is viewed as proportional to both the volume of P3 at the beginning of the current period 

and the size of the resource pool at the beginning of the current period. In addition, since 

plant P3 is able to draw mass from the inaccessible resource pool, plant P3’s growth is 

also proportional to the size of the inaccessible resource pool as well. The growth in P3 

is diminished by the amount of the plant eaten by the non-domesticated herbivores H3 

during the period. Net growth of P3 is shown by (2-107): 

P3 = P3t·(gP3·RPt +gIRPP3·IRPt) –mP3·P3t –H3t·gP3H3·P3t. (2-107) 

The volume of P3 at the end of the current period is equal to the volume at the beginning 

of the current period plus the net natural growth of P3, minus the amount of the plant 

eaten by the non-domesticated herbivores H3 during the period: 

P3t+1 = P3t + P3t·(gP3·RPt +gIRPP3·IRPt) –mP3·P3t –H3t·gP3H3·P3t. (2-108) 

Growth of Non-domesticated Herbivores, H3 

The production (natural growth) of H3 during the current period is viewed as 

proportional to both the amount of H3 at the beginning of the current period and the 

amounts of plants P2 and P3 at the beginning of the current period. In addition, the 

growth of H3 is diminished by the amount of H3 eaten by C2. Net growth is given by 

(109): 

H3 = H3t·(gP2H3·P2t + gP3H3·P3t) –mH3·H3t –C2t·gH3C2·H3t (2-109) 

and the end-of-period stock of H3 is given by: 

H3t+1 = H3t + H3t·(gP2H3·P2t + gP3H3·P3t) –mH3·H3t –C2t·gH3C2·H3t (2-110) 

Growth of Non-domesticated Carnivores, C2 

The mass of C2 increases during the period by the amount of mass that C2 eats of H2 and 

H3; its mass decreases by its natural mortality, mC2·C2t. Net growth of C2 during the 

period is represented by (2-111) and the end-of-period stock of C1 is given by (2-112): 

C2 = C2t·(gH2C2·H2t +gH3C2·H3t) –mC2·C2t (2-111) 

and the end-of-period stock of C2 is given by: 

C2t+1 = C2t + C2t·(gH2C2·H2t +gH3C2·H3t) –mC2·C2t (2-112) 

30
Growth in Human Population and Growth in Human Mass 

Human population at the end of the current period is equal to the human population at the 

beginning of the current period, plus the number of people born during the current period, 

minus human mortality during the period. The human birthrate, η, during the period is 

assumed to be a negative function of the real wage prevailing during the current period. 

The rationale for this assumption is that the real wage represents the opportunity cost of 

opting to remain outside the labor force (at least part-time) for the purpose of rearing 

children. 

Nt+1 = Nt + [η(W/P*) –mN] ·Nt +P1HH +H1HH (2-113) 

– 

where P* denotes a weighted average of the prices the households pay for P1, H1, and IS: 

P*

 = [pP1·P1HH +pH1·H1HH +pIS·ISHH]/ [P1HH +H1HH +ISHH] (2-114) 

The net addition to human mass is equal to human consumption of plants P1 and 

herbivores H1 during the current period, minus human mortality in terms of mass during 

the current period. The net growth in human mass, N^, is given by (2-115): 

N^

t+1 = N^t + P1HH +H1HH –mN·N^t (2-115) 

Mass per human at time t, μt , is given by: 

μt = N^

t/ Nt (2-116) 

Resource Pool 

The size of the resource pool at the end of the current period is equal to its size at the 

beginning of the current period, plus the mass flowing to the resource pool generated by 

plant and animal mortality, minus the amount by which the resource pool declines during 

the current period in direct response to the growth in the domesticated and nondomesticated plants during the current period. In addition the resource pool increases 

over time because of mass transfer from the inaccessible resource pool. This is assumed 

to be in proportion to the size of the inaccessible resource pool. Finally, mass is 

transferred from the accessible resource pool to the inaccessible resource pool in direct 

proportion to the level of production by the IS sector. 

RPt+1 = RPt +mP1·P1t +mP2·P2t +mH1·H1t +mH2·H2t +mC1·C1t +mN·Nt (2-117) 

– P1t·gP1·RPt – P2t·gP2·RPt +rIRPRP·IRPt – λ·IS. 

Inaccessible Resource Pool 

31
The size of the inaccessible resource pool at the end of the current period is equal to its 

size at the beginning of the current period, plus the volume of waste contributed to the 

pool due to the scale of operations of the IS industry, minus the amount that the nondomesticated plants P2 are able to recycle from the inaccessible resource pool during the 

current period and minus the amount of natural decay that moves mass directly from the 

inaccessible resource pool to the resource pool. 

IRPt+1 = IRPt +(θ+λ)·ISHH –P2t·gIRPP2·IRPt – rIRPRP·IRPt. (2-118) 

32
Chapter 3 

Model Solution and Operational Equations 

This chapter provides a model solution based on the assumptions that firms attempt to 

maximize the difference between the sale of their products less the costs (of materials and 

labor), while households maximize their utility (incoming flow of goods balanced against 

leisure). The details and operational equations for the solution are described here. 

Model Solution 

Each firm (P1, H1, IS) seeks to maximize its profit (income minus wages paid and 

purchases of supplies). These industries set their prices in accordance with their current 

forecasts of demand for their products before the market clears. A price setting 

mechanism is used so that trades occur even if markets fail to clear (this is in contrast to a 

general equilibrium setting). Excess inventory is carried to the next time step, and 

shortfalls made up with current inventory. It is assumed that the industrial sector sets the 

money wage rate at the beginning of the current period. In this solution, one of the 

industries, say the IS industry, becomes the dominant industry in the labor market. 

Before any other price is announced, it sets the wage rate, W, at the rate at which it 

anticipates that the amount of labor supplied, net of the sum of the amounts of labor 

demanded by the P1 and H1 industries, will be sufficient to provide the level of labor 

services the IS industry demands. There are several approaches of constructing a solution 

to the model. For example, one could use P1 or H1 as the dominant industries for setting 

wages. However, IS was chosen as the dominant industry because in most modern 

societies, wages tend to be set by the industrial or service sectors rather than the 

agricultural sectors which PI and H1 represent. 

Once the IS industry has announced the wage rate, all industries and the households take 

that wage rate as given. Given the wage rate, the P1, H1 and IS industries decide the 

amount of labor that they want to hire, and the households decide the amount of labor 

they wish to supply to the labor market. It is assumed that the supply of labor exceeds 

the market demand for labor. Thus the amount of labor hired by each industry is the 

amount that each industry demands, given the wage rate announced by the IS industry. 

The three industries simultaneously decide how much labor to hire, decide how much 

output to produce and announce the prices of their respective products, in light of the 

money wage set by the IS industry. At the same time, the IS and H1 industries decide 

their demands for P1. The wage incomes the households receive from the three 

industries, plus their non-wage incomes from the three industries, plus the transfers the 

households receive from the government, plus the revenue from the industries’ fixed 

costs, equals the nominal GDP of the aggregate economy. From (2-95), in real terms, 

GDP and real disposable income, also amounts to the sum: P1HH + H1HH +ISHH. 

Given their real disposable incomes, and the prices the industries set for their products, 

the households then decide their demands for the three products P1, H1, and IS. In each 

33
case, the end of period stocks of P1, H1 and IS correspond to the amount of the product 

produced during the period minus the amount of the item purchased by the various 

sectors. 

Operational Equations 

Linear forms of demand and supply functions are used to simplify calculation. The 

following are consistent with the optimization goal of the firms and the households. 

Firms attempt to maximize the difference between the sale of their products less the costs 

(of materials and labor), while households maximize their utility (incoming flow of goods 

balanced against leisure). The planning horizon in the current model is one time step. 

The model has 19 state variables and 33 outputs. The state variables are given in Table 

3-1. These are in terms of mass, with the exception of the human population. The 

deficits refer to the accumulated difference between the demand for a good and the actual 

quantity of that good delivered in a given time step. Deficits are restricted to being nonpositive. Any surplus of a good increases the inventory of that good. The compartment 

masses are restricted to being non-negative. In addition to the system state variables, 14 

other internal system flows are output. These additional outputs are given in Table 3-2. 

All system flows are restricted to being non-negative. 

Table 3-1 State variables, on a mass basis unless otherwise indicated. 

1 P1 8 C2 15 P1 deficit for HH 

( def P1 HH ) 

2 P2 9 Humans (HH) 16 H1 deficit for HH 

( def H1 HH ) 

3 P3 10 IS 17 IS deficit for HH 

( ISHH def ) 

4 H1 11 RP 18 Number of Humans 

( N HH ) 

5 H2 12 IRP 19 Human per cap mass 

( μ ) 

6 H3 13 P1 deficit for H1 

( def P 1 H1 ) 

7 C1 14 P1 deficit for IS 

( P1IS def ) 

In what follows, the sum of total deficits of either P1, H1 or IS and their current 

inventories represents what the industry must make up in production for past overdemand (or under-supply), in addition to what must be produced to meet the current 

demand. These total deficits to be made up are designated with a star superscript in the 

following ( P1* 

t , H1* t , ISt 

* ): 

34 
P1* 

t = P1t + P1H1t

def + P1IStdef + P1HHtdef 

H1* 

t = H1t + H1HHt

def 

ISt 

* = ISt + ISHHtdef 

Table 3-2. Additional model outputs 

20 RPP1 25 RPIS 30 pP1 

21 P1H2 26 ISHH 31 pH1 

22 P1IS 27 P2H1 32 pIS 

23 P1H1 28 H1C1 33 W 

24 P1HH 29 H1HH 

At a given time step t the IS industry first sets the wage (Wt): 

Wt = aw − cw (ISt 

* − IS) /(λ +θ ) − d w N HH t (3-1) 

In the above, θ and λ convert IS mass to unit quantities. In what follows, conversion 

factors for the other industries are assumed to be unity. The term (⋅) is the end-of-period 

inventory target. Deficits are made up in subsequent time steps, although customers 

already have paid for the goods. 

Based on the wage, the P1, H1 and IS firms set their prices, pt 

(⋅) , and production targets, 

(⋅)t

p , (how much they should produce) in order to maximize profit, given their forecasts 

of demands. All other symbols in the equations below are constant parameters. 

pt

P1 = aP1 + bP1

Wt − cP1 (P1* 

t − P1) (3-2) 

P1t

p = 3d − 3eWt − cP1 p (P1* 

t − P1) (3-3) 

pt

H1 = aH1 + bH1

Wt − cH1 (H1* 

t − H1) (3-4) 

H1t

p = d − eWt − cH1p (H1* 

t − H1) (3-5) 

pt

IS = aIS + bISWt − cIS (ISt 

* − IS ) /(λ + θ ) (3-6) 

ISt

p = d − eWt − cISp (ISt 

* − IS) /(λ +θ ) (3-7) 

Next, the industry and household demands for input goods are calculated. These are 

expressed as flows between compartments. P1H1t is the flow of P1 product to the H1 

industry at time t in units of units, for example: 

35
P1H1t

de = d − eWt − fptP1 − g P1H1 (H1* 

t − H1) (3-8) 

P2H1t

de = k

ˆ (3-9) 

de P1 H1 IS de de de P1HHt = d − kpt + 1

2 mp + 12 np + z(P1HH + H1HH + ISHH ) (3-10) t t t t t 

H1HHt

de = d + 1

2 kpt

P1 − mpH1 + 1

2 npIS + z(P1HH de + H1HH de + ISHH de ) (3-11) t t t t t 

H1 IS de de de ISHHtde = d + 1

2 kp P1 + 12 mp − np + z(P1HH + H1HH + ISHH ) (3-12) t t t t t t 

The last three equations are in terms of units per capita and must be solved 

simultaneously. ISHH de represents the demand for IS by the humans. The flow of goods 

ISHH passes through the human compartment (i.e., these goods are used by humans), but 

does not contribute to the mass of humans. Instead, the mass flows into the inaccessible 

resource pool (IRP). 

Next, the flows that involve labor are calculated. Labor effort is expended to meet 

production targets P1t

p and H1tp as discussed above. The labor goes into fences to keep 

the wild herbivores (H2) and carnivores (C1) from eating stocks: 

P1H 2t = g RPP1P1t RPt − P1t

p − mP1P1t (3-13) 

H1C1t = P1H1t + P2H1t − H1t

p − mH1H1t (3-14) 

This is completed by calculating IS industry demands flows. IS demand and supply is in 

individual units and thus requires a conversion factor to mass for each input. This is true 

of all goods, but conversions for P1 and H1 are assumed to be unity (one unit requires 1 

unit of mass). 

P1ISt

de = IStp

θ (3-15) 

RPISt

de = IStp λ (3-16) 

For the next time step, in the following, gi are growth rates and mi are mortality rates for 

the individual species i. All state variables are in terms of mass, except for NHH t . Checks 

are done so that all flows are positive, and that compartments maintain a non-negative 

mass (except the deficit state variables). This latter requires that flows sometimes do not 

meet demands. Such deficits are accumulated in the deficit state variables and contribute 

to the production functions as outlined above. These deficit state variables are given by: 

P1H1def = P1H1def + P1H1 − P1H1de (3-17) t+1 t t t 

def def de P1ISt+1 = P1ISt + P1ISt − P1ISt (3-18) 

def def de P1HHt +1 = P1HHt + P1HHt − P1HHt N HH t (3-19) 

36
def def de H1HHt +1 = H1HHt + H1HHt − H1HHt N HH t (3-20) 

ISHH def = ISHH def + ISHH − (θ + λ)ISHH de * N (3-21) t+1 t t t HH t 

The above simply track the difference between what was actually delivered versus 

demanded. For example, P1IS refers to the flow from P1 to IS, which may be less or 

greater than P1ISt

de , depending on the situation at that time step. 

Excess surplus is carried over to the next time step, or is used to make up deficits flows 

from previous time steps. The remaining state equations are given as follows. 

P1t+1 = P1t + P1t (g RPP1 RPt − mP1 ) − P1H1t − P1H 2t − P1HHt − P1ISt (3-22) 

100 P2t+1 = P2t + P2t (g RPP2 RPt + rIRPP2 2 IRPt − mP2 − g P2H 2 H 2t − g P2H 3 H 3t ) − P2H1t 100 + IRP 

(3-23) 

100 P3t+1 = P3t + P3t (g RPP3 RPt + rIRPP3 IRPt − mP3 − g P3H 3 H 3t ) (3-24) 100 + IRP2 

H1t +1 = H1t + P1H1t + P2H1t − mH1 H1t − H1C1t − H1HHt (3-25) 

H 2t+1 = H 2t + P1H 2t + H 2t (g P2H 2 P2t − mH 2 − g H 2C1C1t − g H 2C 2C2t ) (3-26) 

H 3t+1 = H 3t + H 3t (g P2H 3 P2t + g P3H 3 P3t − mH 3 − g H 3C 2C2t ) (3-27) 

C1t+1 = C1t + H1C1t + C1t (g H 2C1H 2t − mC1 ) (3-28) 

C2t+1 = C2t + C2t (g H 2C 2 H 2t + g H 3C 2 H 3t − mC 2 ) (3-29) 

HHt+1 = HHt + P1HHt + H1HHt − mHH N HH t μt (3-30) 

N = N + (η(

Wt 

Pt 

* ) − mHH −φ(μt − μ I )

2

)N HH t (3-31) HH t +1 HH t 

μt+1 = HHt+1 

N HH t+1 

(3-32) 

ISt+1 = ISt + P1ISt + RPISt − ISHHt (3-33) 

37
RPt+1 = RPt + mP1 P1t + mP2 P2t + mP3 P3t + mH1 H1t + mH 2 H 2t + mH 3 H 3t 

+ mC1C1t + mC 2C2t + mHH HHt + mIRPRP IRPt − RPt (g RPIRP IRPt (3-34) 

+ g RPP1 P1t + g RPP2 P2t + g RPP3 P3t ) − RPISt 

IRPt +1 = IRPt + ISHHt − IRPt (rIRPP2 P2t + rIRPP3 P3t + mIRPRP ) (3-35) 

Where 

Pt 

* = ( ptP1 P1HHt + ptH1 H1HHt + ptIS ISHHt ) /(P1HHt + H1HHt + ISHHt ) (3-36) 

η(

Wt 

P* ) = aη − bη (

Wt 

P* ) 

1

2 (3-37) t t 

All parameters and their nominal values are given in Table 3-3 along with the 

corresponding state variable initial conditions in Table 3-4. 

38
Table 3-3. Parameters and their nominal values

* 

Ecological Parameters Economic Parameters 

1. gRPP2=2.861325e-2 1. Aw=4.832249e-1 32. mP1HH=1.997642e-4 

2. gP2H2=2.934352e-2 2. Cw=1.357181e-1 33. nP1HH=7.776176e-5 

3. gP2H3=0.042 3. aP1=1.0 34. dH1HH=1.9108e-004 

4. gRPP3=5.732733e-3 4. bP1=1.0 35. zH1HH=1.4745e-001 

5. gP3H3=3.131235e-1 5. cP1=7.737088e-1 36. kH1HH=2.1204e-004 

6. gH2C1=9.174907e-1 6. aP1p=5.732311e-1 37. mH1HH=3.9953e-004 

7. gH2C2=1.312728e-1 7. bP1p=1.497375e-1 38. nH1HH=7.7762e-005 

8. gH3C2=2.938371e-1 8. cP1p=3.380538e-2 39. dISHH=1.9108e-004 

9. rIRPP2=1.081656e-2 9. aH1=7.524099e-1 40. zISHH=1.4745e-001 

10. rIRPP3=0.9 10. bH1=0.001 41. kISHH=2.1204e-004 

11. mP2=4.932829e-1 11. cH1=2.527165e-1 42. mISHH=1.9976e-004 

12. mP3=4.658138e-1 12. aH1p=1.910770e-1 43. nISHH=1.5552e-004 

13. mH2=0.001 13. bH1p=4.991250e-2 44. K^=3.0000e-001 

14. mH3=4.903092e-1 14. cH1p=9.926235e-1 45. θ (Theta)=1.0199e-001 

15. mC1=2.302639e-1 15. aIS=6.081040e-1 46. λ (Lambda)=6.7668e-001 

16. mC2=4.286472e-1 16. bIS=2.972103e-1 47. gRPP1=0.09 

20. gP1H2=0.1 17. cIS=0.001 48. mP1=1.018295e-3 

21. gH1C1=0.2 18. aISp=1.910770e-1 49. mH1=9.838862e-3 

19. bISp=4.991250e-2 50. mHH=0.22 

20. cISp=5.646818e-1 51. P 1 (P1bar)mass=0 

21. dP1H1=1.910770e-4 52. H 1 (H1bar)mass=0.4 

22. eP1H1=4.991250e-2 53. I S (ISbar)mass=0 

23. fP1H1=4.240783e-1 54. Dw=4.507354e-6 

24. gP1H1=1.9 55. η (Eta)a=1.5000e+000 

29. dP1HH=1.910770e-4 56. η (Eta)b=8.3333e-001 

30. zP1HH=1.474467e-1 57. η (Eta)c=1.6667e-002 

31. kP1HH=4.240783e-4 58. n (Phi) =1.0000e+001 

59. Idealpercapmass=4.0556e003 

*

 Note, these are grouped by their context in the software, not necessarily by number. 

39 
Table 3-4. Corresponding state variable initial conditions. 

1 P1 1.55031758212386 8 C2 1.32270575776743 15 def P1 HH 0 

2 P2 10.00000000000000 9 HH 0.45073543305286 16 def H1 HH 0 

3 P3 1.46577689115230 10 IS 0.10089700803626 17 ISHH def 0 

4 H1 0.14002360565832 11 RP 19.14708496882624 18 N HH 10 

5 H2 0.25097142084356 12 IRP 0.86413620441690 19 μ 0.004055586 

60633 

6 H3 1.34666957783805 13 def P 1 H1 0 

7 C1 0.12948354945840 14 P1IS def 0 

40
Chapter 4 

Summary 

Sustainability is fundamentally an effort to create and maintain a regime in which the 

human population and its necessary energy and material consumption can be supported 

indefinitely by the biological system of the Earth. The model system presented in this 

report consists of a foodweb with an integrated industrial sector, agricultural production, 

a very simple economy, and the rudiments of a legal system. As such, it represents a 

promising first step in understanding the relationship between ecology, economics, 

technology, and society. The purpose of this report is to document the model so as to 

make it accessible to the larger user community in the hopes of stimulating further 

research work. 

The preliminary model described here was developed in stages. First, a basic ecological 

model was created with the total mass distributed among all compartments within a 

foodweb. The compartments were then differentiated in terms of human control over 

them (by identifying species as being either domesticated or non-domesticated, assigning 

property rights to the mass in certain compartments, and indicating intentional changes to 

the natural flows of mass, such as fences). Finally, a generic industrial process (with 

private property rights) was added that diverts resources to human (non-food) 

consumption and use. 

From an ecological perspective, the model contains two resource pools (one biologically 

accessible and one inaccessible), three primary producers (P1, P2, and P3); three 

herbivores (H1, H2, and H3), two carnivores (C1 and C2), and humans (HH). The system 

flows are specified in terms of mass. The system is closed to mass (i.e., mass is 

conserved) and open to a non-limiting source of energy. 

From an economic perspective the model contains human households, an industrial sector 

(IS), and two private agricultural firms: one a producer of plants (P1) and one a producer 

of herbivores (H1). The households are the ultimate owners of the factors of production 

used to produce the goods that are traded in explicit markets. In addition to markets for 

the three goods (P1, H1, IS), there is a labor market. 

The optimal economic behavior of industry, government, and households is as follows. 

The P1 industry applies a variable amount of labor and a fixed amount of capital to 

transform the mass of P1 that would otherwise grow naturally into a marketable product. 

The P1 industry hires labor to reduce the consumption of P1 by H2. The H1 industry 

buys P1, labor and grazing rights to P2 and sells output to the households. The H1 

industry hires labor to limit the amount of H1 eaten by C1. The IS industry buys P1 and 

combines it with mass from the resource pool (RP) using variable labor and a fixed 

amount of capital to produce IS, which it sells to the households. The government sector 

receives revenue from the grazing rights it grants to the H1 industry and transfers this 

revenue to the households. The objective of the household sector is to maximize current 

41
utility, which is assumed to be a positive increasing function of consumption per capita of 

P1, H1, IS and leisure. This economic activity is reflected in the nominal gross national 

product (GDP), which is the total market value of all units of final goods produced in the 

economy during the period, which also corresponds to the sum of all spending on final 

goods plus the values of the changes in inventories of all goods. 

The production (natural growth) of non-domesticated species during the period is 

governed by biological rules, but is influenced by human economic behavior (i.e., fences 

to protect domesticated species). 

The human population at the end of the current period is equal to the human population at 

the beginning of the current period, plus the number of people born during the current 

period, minus human mortality during the period. The human birthrate during the period 

is assumed to be a negative function of the real wage prevailing during the current period. 

The net addition to human mass is equal to human consumption of plants P1 and 

herbivores H1 during the current period, minus human mortality in terms of mass during 

the current period. 

The size of the resource pool at the end of the current period is equal to its size at the 

beginning of the current period, plus the mass flowing to the resource pool generated by 

plant and animal mortality, minus the amount by which the resource pool declines during 

the current period in direct response to the growth in the domesticated and nondomesticated plants during the current period. In addition the resource pool increases 

over time because of mass transfer from the inaccessible resource pool. Mass is 

transferred from the resource pool to the inaccessible resource pool in direct proportion to 

the level of production by the IS sector. 

The size of the inaccessible resource pool at the end of the current period is equal to its 

size at the beginning of the current period, plus the volume of waste contributed to the 

pool due to the scale of operations of the IS industry, minus the amount that the nondomesticated plants P1 and P2 are able to recycle from the inaccessible resource pool 

during the current period and minus the amount of natural decay that moves mass directly 

from the inaccessible resource pool to the resource pool. 

This study is a preliminary one and did not include all possible parameters or model behaviors. For example, it is known that changes in the recycling parameter for P2 and P3 

can cause the system to crash entirely. The solution used in this report consisted of the 

following iterative steps. 

• The industrial sector sets the wage rate. The assumption is that this sector 

dominates the labor market. 

• Based on the wage, all industries set their prices and their production targets 

according to their own utility functions and internal models of demand of their 

products. 

• Humans determine their demands for goods. 

• Industries determine their demands for goods and labor. 

42
• Checks are done for internal consistencies of flows (to be sure they meet 

positivity constraints on flows and compartment masses). 

• The next step is taken (flows are transferred) for both the economic and 

ecological parts of the model. 

Some preliminary simulations have been conducted using this model (Cabezas et al. 

2006). However, the modle behavior has not been explored in detail. Future uture work 

entails further refinement of the model and methods of analysis. Several options are 

being considered for expanding the model and attempting to use it to suggest or evaluate 

possible policy changes that would allow all species to continue while maintaining a 

viable economic environment for the human population. Some of these options include 

adding economic and regulatory tools such a fee on waste generation, futher legal 

protections of wild species, and inclusion of environmental quality valuation. The overall 

general goal is to use the model as rough generic screening tool for sustainable 

environmental management strategies. This might include the use of economic 

incentives and regulation. 

43
Appendix A 

Symbols Used in the Report

C1t = the stock of C1 at the beginning of the period 

cP1H2 = consumption of P1 per unit of H2 (depends on the restrictions imposed by P1 

industry) 

gi = growth parameter for species i 

gP1 = growth parameter of P1 

H1s

 = H1 industry’s supply of H1 during the current period 

H1t = stock of H1 at the beginning of the period 

H 1 = the H1 industry’s desired end of period inventory of H1 

H1HH = the amount of H1 consumed by the households 

(H1HH)d

 = household sector’s demand for H1 

H2t = stock of H2 at the beginning of the period 

hH1 = average number of hours employees work in the H1 industry producing H1 

h*

H1 = average number of hours employees are used in the H1 industry to reduce the 

amount of H1 eaten by C1 

hIS = average number of hours employees work in the IS industry producing IS 

hNd = labor hours demanded 

hP1 = average number of hours employees work in the P1 industry producing P1 

h*

P1 = average number of hours employees are used in the P1 industry to reduce the 

amount of P1 eaten by H2 

hs

 = average number of hours the households would like to work during the period (for all 

employers combined) 

ISt = the initial stock of IS (in units of goods) held by the IS industry 

(θ+λ)·ISt = the initial stock of IS (in units of mass) held by the IS industry 

44
I S = the ending inventory of IS (in units of goods) that the IS industry plans for the end 

of the period 

(ISHH)s

 = the amount of IS the IS industry desires to supply to the households 

(ISHH)d

 = the household sector’s demand for IS 

KH1

t = stock of physical capital held by the H1 industry at the beginning of the period 

KIS

t = stock of physical capital held by the IS industry at the beginning of the period 

KP1

t = stock of physical capital held by the P1 industry at the beginning of the period 

K^ = P2H1 (the amount of P2 the government permits the H1 industry to take) 

mi = mortality parameter for species i 

mP1 = mortality parameter of P1 

mH1 = mortality parameter of H1 

NS

 = net supply of labor function 

Nt = human population at time t 

N^ 

t = total human mass at time t 

P*

 = (weighted) average price of consumer goods 

P1H1 = the amount of P1 consumed by the H1 industry 

(P1H1)d

 = the H1 industry’s demand for P1 

P1HH = the amount of P1 consumed by the households 

(P1HH)d

 = the household sector’s total demand for P1 

(P1IS)d

 = the IS industry’s demand for P1 

P1s = P1 industry’s supply of P1 during the period 

P1t = stock of P1 at the beginning of the period 

P 1 = desired inventory of P1 for the end of the period 

45
P2H1 = amount of P2 the government permits the H1 industry to take 

pH1 = price of H1 announced by the H1 industry 

pP1 = price of P1 announced by the P1 industry 

RPt = mass of the resource pool at the beginning of the period 

Tr = the transfer payments received by the households from the government sector 

W = money wage rate 

WhNd = total wage bill (wage rate times number of labor hours) 

η(•) = human birthrate 

μt = average human mass per capita at time t = N^ 

t / N t 

Π = total dividend income received by the households from all industries 

Π*

P1 = profits of the P1 industry 

 = dividends plus additions to retained earnings for P1 industry 

 = dividends plus unintended real investment in inventories in the P1 industry 

Π** = total non-wage income earned by all industries 

 = profits plus rental income to fixed factors earned by all industries 

n = weighting factor for the difference between the current per capita weight of a human 

and the ideal weight 

λ = the amount of RP necessary to produce a unit of IS 

θ = the amount of P1 necessary to produce a unit of IS 

θ·IS = the amount of P1 used by the IS industry 

(θ+λ)·IS = the amount of mass appearing in the IS produced during the period 

(θ+λ)·ISHH = the amount of mass the households buy during the current period in the 

form of IS 

ωH1 = a vector of parameters other than pP1 that the P1 industry believes affects the H1 

industry’s demand for P1 

ωIS = a vector of parameters other than pP1 that the P1 industry believes affects the IS 

industry’s demand for P1 

46
ωISHH = a vector of parameters other than pIS that the IS industry believes affects the 

households’ demand for IS 

ωHH = a vector of parameters other than pP1 that the P1 industry believes affects the 

households’ demand for P1 

ωH1HH = a vector of parameters other than pH1 that the H1 industry believes affects the 

households’ demand for H1 

47
Run price_setting_initialize10.m to initialize parameters 

price_settingSv8 x

S-Function To Workspace1 

t 

Clock To Workspace 

Appendix B 

SIMULINK Graphical Model 

Figure B-1. Root level SIMULINK model (file price_setting2.mdl as opened in SIMULINK) 

The S-Function contains a pointer to the file price_settingSv8.m 

It is a masked system, with the mask defining the parameter names (values loaded via 

initialization file price_setting_initialization10.m). 

To Workspace1 defines a MATLAB workspace variable x which holds all the model 

output. 

To Workspace defines a MATLAB workspace variable t which holds the time variable. 

48
Appendix C 

SIMULINK Code 

Text contained in and defining the graphical model of Appendix 1. Filename: 

price_setting2.mdl 

Model {

Name "price_setting2"

Version 6.1 

MdlSubVersion 0 

GraphicalInterface {

NumRootInports 0 

NumRootOutports 0 

ParameterArgumentNames "" 

ComputedModelVersion "1.151" 

NumModelReferences 0 

NumTestPointedSignals 0 

}

SavedCharacterEncoding "ibm-5348_P100-1997" 

SaveDefaultBlockParams on 

SampleTimeColors off 

LibraryLinkDisplay "none" 

WideLines off 

ShowLineDimensions off 

ShowPortDataTypes off 

ShowLoopsOnError on 

IgnoreBidirectionalLines off

ShowStorageClass off 

ShowTestPointIcons on 

ShowViewerIcons on 

SortedOrder off 

ExecutionContextIcon off 

ShowLinearizationAnnotations on 

RecordCoverage off 

CovPath "/"

CovSaveName "covdata" 

CovMetricSettings "dw" 

CovNameIncrementing off 

CovHtmlReporting on 

covSaveCumulativeToWorkspaceVar on

CovSaveSingleToWorkspaceVar on

CovCumulativeVarName "covCumulativeData" 

CovCumulativeReport off 

CovReportOnPause on 

ScopeRefreshTime 0.035000 

OverrideScopeRefreshTime on

DisableAllScopes off 

DataTypeOverride "UseLocalSettings"

MinMaxOverflowLogging "UseLocalSettings"

MinMaxOverflowArchiveMode "Overwrite" 

BlockNameDataTip off 

BlockParametersDataTip off 

49
 BlockDescriptionStringDataTip off 

ToolBar on 

StatusBar on 

BrowserShowLibraryLinks off

BrowserLookUnderMasks off 

Created "Tue Aug 12 09:54:50 2003"

UpdateHistory "UpdateHistoryNever"

ModifiedByFormat "%<Auto>" 

LastModifiedBy "CPAWLOWS" 

ModifiedDateFormat "%<Auto>" 

LastModifiedDate "Mon Mar 28 14:46:08 2005" 

ModelVersionFormat "1.%<AutoIncrement:151>" 

ConfigurationManager "None" 

LinearizationMsg "none" 

Profile off 

ParamWorkspaceSource "MATLABWorkspace"

AccelSystemTargetFile "accel.tlc" 

AccelTemplateMakefile "accel_default_tmf" 

AccelMakeCommand "make_rtw" 

TryForcingSFcnDF off 

ExtModeBatchMode off 

ExtModeEnableFloating on 

ExtModeTrigType "manual" 

ExtModeTrigMode "normal" 

ExtModeTrigPort "1" 

ExtModeTrigElement "any"

ExtModeTrigDuration 1000 

ExtModeTrigDurationFloating "auto"

ExtModeTrigHoldOff 0 

ExtModeTrigDelay 0 

ExtModeTrigDirection "rising"

ExtModeTrigLevel 0 

ExtModeArchiveMode "off" 

ExtModeAutoIncOneShot off 

ExtModeIncDirWhenArm off 

ExtModeAddSuffixToVar off 

ExtModeWriteAllDataToWs off 

ExtModeArmWhenConnect on 

ExtModeSkipDownloadWhenConnect off

ExtModeLogAll on 

ExtModeAutoUpdateStatusClock on

BufferReuse on 

StrictBusMsg "None" 

ProdHWDeviceType "32-bit Generic" 

ShowModelReferenceBlockVersion off 

ShowModelReferenceBlockIO off 

Array {

Type "Handle" 

Dimension 1 

Simulink.ConfigSet {

$ObjectID 1 

Version "1.0.4" 

Array {

Type "Handle" 

Dimension 7 

 Simulink.SolverCC { 

$ObjectID 2 

50
 Version "1.0.4" 

StartTime "0.0" 

StopTime "200" 

AbsTol "auto" 

FixedStep "auto" 

InitialStep "auto" 

MaxNumMinSteps "-1" 

MaxOrder 5 

ExtrapolationOrder 4 

NumberNewtonIterations 1 

MaxStep "auto" 

MinStep "auto" 

RelTol "1e-3" 

SolverMode "Auto" 

Solver "FixedStepDiscrete"

SolverName "FixedStepDiscrete"

ZeroCrossControl "UseLocalSettings"

PositivePriorityOrder off 

AutoInsertRateTranBlk off 

SampleTimeConstraint "Unconstrained" 

RateTranMode "Deterministic" 

}

 Simulink.DataIOCC { 

$ObjectID 3 

Version "1.0.4" 

Decimation "1" 

ExternalInput "[t, u]"

FinalStateName "xFinal" 

InitialState "xInitial" 

LimitDataPoints on 

MaxDataPoints "1000" 

LoadExternalInput off 

LoadInitialState off 

SaveFinalState off 

SaveFormat "Array"

SaveOutput on 

SaveState off 

SignalLogging on 

SaveTime on 

StateSaveName "xout" 

TimeSaveName "tout" 

OutputSaveName "yout"

SignalLoggingName "sigsOut"

OutputOption "RefineOutputTimes"

OutputTimes "[]"

Refine "1" 

}

 Simulink.OptimizationCC { 

$ObjectID 4 

Version "1.0.4" 

BlockReduction on 

BooleanDataType on 

ConditionallyExecuteInputs on

ConditionalExecOptimization "on_for_testing"

InlineParams off 

InlineInvariantSignals on 

OptimizeBlockIOStorage on 

51
 BufferReuse on 

EnforceIntegerDowncast

ExpressionFolding

FoldNonRolledExpr

LocalBlockOutputs

ParameterPooling

RollThreshold 

on 

on 

on 

on 

on 

5 

SystemCodeInlineAuto

StateBitsets 

off 

off 

DataBitsets off 

UseTempVars off 

ZeroExternalMemoryAtStartup on

ZeroInternalMemoryAtStartup on

InitFltsAndDblsToZero on 

NoFixptDivByZeroProtection off

OptimizeModelRefInitCode off

LifeSpan "inf" 

}

 Simulink.DebuggingCC { 

$ObjectID 5 

Version "1.0.4" 

RTPrefix "error" 

ConsistencyChecking "none" 

ArrayBoundsChecking "none" 

AlgebraicLoopMsg "warning"

ArtificialAlgebraicLoopMsg "warning"

CheckSSInitialOutputMsg on

CheckExecutionContextPreStartOutputMsg off

CheckExecutionContextRuntimeOutputMsg off 

SignalResolutionControl "TryResolveAllWithWarning"

BlockPriorityViolationMsg "warning"

MinStepSizeMsg "warning"

SolverPrmCheckMsg "none" 

InheritedTsInSrcMsg "warning"

DiscreteInheritContinuousMsg "warning"

MultiTaskDSMMsg "warning"

MultiTaskRateTransMsg "error" 

SingleTaskRateTransMsg "none" 

TasksWithSamePriorityMsg "warning"

CheckMatrixSingularityMsg "none"

IntegerOverflowMsg "warning"

Int32ToFloatConvMsg "warning"

ParameterDowncastMsg "error" 

ParameterOverflowMsg "error" 

ParameterPrecisionLossMsg "warning"

UnderSpecifiedDataTypeMsg "none"

UnnecessaryDatatypeConvMsg "none"

VectorMatrixConversionMsg "none"

InvalidFcnCallConnMsg "error" 

FcnCallInpInsideContextMsg "Use local settings"

SignalLabelMismatchMsg "none" 

UnconnectedInputMsg "warning"

UnconnectedOutputMsg "warning"

UnconnectedLineMsg "warning"

SFcnCompatibilityMsg "none" 

UniqueDataStoreMsg "none" 

RootOutportRequireBusObject "warning" 

52
 AssertControl "UseLocalSettings"

EnableOverflowDetection off 

ModelReferenceIOMsg "none" 

ModelReferenceVersionMismatchMessage "none"

ModelReferenceIOMismatchMessage "none"

ModelReferenceCSMismatchMessage "none"

ModelReferenceSimTargetVerbose off

UnknownTsInhSupMsg "warning"

ModelReferenceDataLoggingMessage "warning"

ModelReferenceSymbolNameMessage "warning"

}

 Simulink.HardwareCC { 

$ObjectID 6 

Version "1.0.4" 

ProdBitPerChar 8 

ProdBitPerShort 16 

ProdBitPerInt 32 

ProdBitPerLong 32 

ProdIntDivRoundTo "Undefined" 

ProdEndianess "Unspecified"

ProdWordSize 32 

ProdShiftRightIntArith on 

ProdHWDeviceType "32-bit Generic" 

TargetBitPerChar 8 

TargetBitPerShort 16 

TargetBitPerInt 32 

TargetBitPerLong 32 

TargetShiftRightIntArith on

TargetIntDivRoundTo "Undefined" 

TargetEndianess "Unspecified"

TargetWordSize 32 

TargetTypeEmulationWarnSuppressLevel 0

TargetPreprocMaxBitsSint 32

TargetPreprocMaxBitsUint 32

TargetHWDeviceType "Specified"

TargetUnknown on 

ProdEqTarget on 

}

 Simulink.ModelReferenceCC { 

$ObjectID 7 

Version "1.0.4" 

UpdateModelReferenceTargets "IfOutOfDateOrStructuralChange"

CheckModelReferenceTargetMessage "error"

ModelReferenceNumInstancesAllowed "Multi" 

ModelReferencePassRootInputsByReference on

ModelReferenceMinAlgLoopOccurrences off

}

 Simulink.RTWCC { 

$BackupClass "Simulink.RTWCC" 

$ObjectID 8 

Version "1.0.4" 

SystemTargetFile "grt.tlc"

GenCodeOnly off 

MakeCommand "make_rtw" 

TemplateMakefile "grt_default_tmf"

GenerateReport off 

SaveLog off 

53
 RTWVerbose on 

RetainRTWFile off 

ProfileTLC off 

TLCDebug

TLCCoverage

TLCAssert 

off 

off 

off 

ProcessScriptMode

ConfigurationMode

ConfigAtBuild 

"Default" 

"Optimized"

off 

IncludeHyperlinkInReport off

LaunchReport off 

Array {

Type "Handle" 

Dimension 2 

Simulink.CodeAppCC {

$ObjectID 9 

Version "1.0.4" 

ForceParamTrailComments off 

GenerateComments on 

IgnoreCustomStorageClasses on

IncHierarchyInIds off 

MaxIdLength 31 

PreserveName off 

PreserveNameWithParent off 

ShowEliminatedStatement off 

IncAutoGenComments off 

SimulinkDataObjDesc off 

SFDataObjDesc off 

IncDataTypeInIds off 

PrefixModelToSubsysFcnNames on

CustomSymbolStr "$R$N$M" 

MangleLength 1 

DefineNamingRule "None" 

ParamNamingRule "None" 

SignalNamingRule "None" 

InsertBlockDesc off 

SimulinkBlockComments on 

EnableCustomComments off 

InlinedPrmAccess "Literals" 

ReqsInCode off 

}

Simulink.GRTTargetCC {

$BackupClass "Simulink.TargetCC"

$ObjectID 10 

Version "1.0.4" 

TargetFcnLib "ansi_tfl_tmw.mat" 

GenFloatMathFcnCalls "ANSI_C" 

UtilityFuncGeneration "Auto" 

GenerateFullHeader on 

GenerateSampleERTMain off 

IsPILTarget off 

ModelReferenceCompliant off

IncludeMdlTerminateFcn on 

CombineOutputUpdateFcns off

SuppressErrorStatus off 

IncludeFileDelimiter "Auto" 

ERTCustomFileBanners off 

54
 SupportAbsoluteTime on 

LogVarNameModifier "rt_" 

MatFileLogging off 

MultiInstanceERTCode off 

SupportNonFinite on 

SupportComplex on 

PurelyIntegerCode off 

SupportContinuousTime on 

SupportNonInlinedSFcns on 

ExtMode off 

ExtModeStaticAlloc off 

ExtModeTesting off 

ExtModeStaticAllocSize 1000000 

ExtModeTransport 0 

ExtModeMexFile "ext_comm" 

RTWCAPISignals off 

RTWCAPIParams off 

RTWCAPIStates off 

GenerateASAP2 off 

}

PropName "Components"

}

}

PropName "Components"

}

Name "Configuration"

SimulationMode "normal" 

CurrentDlgPage "Solver" 

}

PropName "ConfigurationSets"

}

Simulink.ConfigSet {

$PropName "ActiveConfigurationSet"

$ObjectID 1 

}

BlockDefaults {

Orientation "right"

ForegroundColor "black" 

BackgroundColor "white" 

DropShadow off 

NamePlacement "normal" 

FontName "Helvetica" 

FontSize 10 

FontWeight "normal" 

FontAngle "normal" 

ShowName on 

}

BlockParameterDefaults {

Block {

BlockType Clock 

DisplayTime off 

}

Block {

BlockType "S-Function" 

FunctionName "system"

SFunctionModules "''" 

PortCounts "[]" 

55
 }

Block {

BlockType

VariableName 

MaxDataPoints 

ToWorkspace

"simulink_output"

"1000" 

Decimation "1" 

SampleTime

FixptAsFi

}

}

AnnotationDefaults {

HorizontalAlignment

VerticalAlignment

ForegroundColor

BackgroundColor

DropShadow

FontName 

"0" 

off 

"center" 

"middle" 

"black" 

"white" 

off 

"Helvetica" 

FontSize 10 

FontWeight

FontAngle

}

LineDefaults {

FontName 

"normal" 

"normal" 

"Helvetica" 

FontSize 9 

FontWeight

FontAngle

} 

"normal" 

"normal" 

System {

Name "price_setting2"

Location [368, 122, 948, 423]

Open on 

ModelBrowserVisibility off 

ModelBrowserWidth 200 

ScreenColor "white" 

PaperOrientation "landscape"

PaperPositionMode "auto" 

PaperType "usletter" 

PaperUnits "inches" 

ZoomFactor "100" 

ReportName "simulink-default.rpt"

Block {

BlockType Clock 

Name "Clock" 

Position [160, 175, 180, 195]

Decimation "10" 

}

Block {

BlockType "S-Function" 

Name "S-Function" 

Ports [0, 1]

Position [85, 80, 235, 110]

FunctionName "price_settingSv8"

Parameters 

"ecolparams,econparams,ic,MH_updown,Lamb_updown,"

"Thet_updown, tlo, thi"

MaskPromptString "Human Mortality|Lambda, RP input for IS

industr" 

56
"y|Theta, P1 input for IS industry|tlo, time to start change|thi, time

to end " 

"change"

MaskStyleString "popup(unchanging|increasing to not more

than 2 " 

"times|decreasing to not less than 1/10th),popup(unchanged|increasing

w/o boun"

"d|decreasing to not less than 1/10th),popup(unchanging|increasing w/o

bound|d"

"ecreasing to not less than 1/10th),edit,edit"

MaskTunableValueString "off,off,off,off,off"

MaskCallbackString "||||"

MaskEnableString "on,on,on,on,on"

MaskVisibilityString "on,on,on,on,on"

MaskToolTipString "on,on,on,on,on"

MaskVarAliasString ",,,,"

MaskVariables 

"MH_updown=@1;Lamb_updown=@2;Thet_updown=@3;tlo="

"@4;thi=@5;"

MaskIconFrame on 

MaskIconOpaque

MaskIconRotate 

on 

"none" 

MaskIconUnits "autoscale" 

MaskValueString

MaskTabNameString

}

Block {

BlockType

Name 

Position 

VariableName 

"unchanging|unchanged|unchanging|10|10"

",,,," 

ToWorkspace

"To Workspace"

[325, 170, 385, 200]

"t" 

MaxDataPoints "inf" 

SampleTime

SaveFormat 

}

Block {

BlockType

Name 

Position 

VariableName 

"1" 

"Array" 

ToWorkspace

"To Workspace1"

[335, 80, 395, 110]

"x" 

MaxDataPoints "inf" 

SampleTime

SaveFormat 

}

Line {

SrcBlock 

"1" 

"Array" 

"Clock" 

SrcPort 1 

DstBlock 

DstPort 

"To Workspace"

1 

}

Line {

SrcBlock "S-Function" 

SrcPort 1 

DstBlock 

DstPort 

"To Workspace1"

1 

}

Annotation { 

57
 Name "Run price_setting_initialize10.m to initialize

" 

"parameters"

Position [258, 26]

}

}

} 

58
Appendix D 

MATLAB Code 

The first file is an initialization file, price_setting_initialize10.m, which initializes the 

parameters and initial state conditions. This is a convoluted way of loading the 

parameters listed above. 

The next file is the core of the model, written as an S function. The file name is 

price_settingSv8.m 

Initialization file (price_setting_initialize10.m). Run before running SIMULINK model 

% price_setting_initialize10.m

%initialization file for price_settingSv8.m

%

 ecolparams=[2.861325308717139e-002 %found 30 Aug 2004

1.734351784507601e-002

 3.265224224356716e-002

 5.732733145333690e-003

 3.131235080180339e-001

 9.174906638699569e-001

 3.327275637705977e-001

 2.350696927025789e-001

 5.081656269272850e-002

 8.688268022462951e-001

 4.932828648894071e-001

 4.658138102587053e-001

 1.000000000000000e-003

 5.030915711147533e-001

 2.302639355221487e-001

 4.286472149788781e-001

 0.0%6.602439763447394e-006 %background feedback of IRP to RP set to

% zero

 1.233437625619434e+000

 1.000000000000000e+001

 .1 

.2]; 

%gP1H2 %when humans die off

%gH1C1 

ecolparams(10)=.9;

ecolparams(7)=1.3127275637705977e-001; 

ecolparams(14)=4.9030915711147533e-001;

ecolparams(3)=3.7e-2;

ecolparams(9)=1.081656269272850e-002;

ecolparams(2)=2.934351784507601e-002;

ecolparams(3)=4.2e-2; %6 Oct 2004 

59
% ecological params =

% [gRPP2;gP2H2;gP2H3;gRPP3;gP3H3;gH2C1;gH2C2;gH3C2;

% rIRPP2;rIRPP3;mP2;mP3;mH2;mH3;mC1;mC2;mIRPRP;RPIRP;P2;gP1H2;gH1C1];

% Here, P2 is a

% choice variable, controlling the total system mass

%

% gRPP2=ecolparams(1);

% 

% gP2H2=ecolparams(2);gP2H3=ecolparams(3);gRPP3=ecolparams(4);

% gP3H3=ecolparams(5);

% gH2C1=ecolparams(6);gH2C2=ecolparams(7);gH3C2=ecolparams(8);

% rIRPP2=ecolparams(9);rIRPP3=ecolparams(10);

% mP2=ecolparams(11);mP3=ecolparams(12);mH2=ecolparams(13);

% mH3=ecolparams(14);mC1=ecolparams(15);mC2=ecolparams(16);

% mIRPRP=ecolparams(17);RPIRP=ecolparams(18); %not used in full econ-

% ecol model

% P2=ecolparams(19);

% gP1H2=ecolparams(20);gH1C1=ecolparams(21); %13 Aug 2004 natural

% predation of P1 and H1

%

% economic parameters, found using ga_find_econV1.m:

%

% aw=econparams(1);cw=econparams(2);

% aP1=econparams(3);bP1=econparams(4);cP1=econparams(5);

% aP1p=econparams(6);bP1p=econparams(7);cP1p=econparams(8);

% aH1=econparams(9);bH1=econparams(10);cH1=econparams(11);

% aH1p=econparams(12);bH1p=econparams(13);cH1p=econparams(14);

% aIS=econparams(15);bIS=econparams(16);cIS=econparams(17);

% aISp=econparams(18);bISp=econparams(19);cISp=econparams(20);

% 

% dP1H1=econparams(21);eP1H1=econparams(22);fP1H1=econparams(23);

% gP1H1=econparams(24);

% 

% *****dP1IS=econparams(25);eP1IS=econparams(26);fP1IS=econparams(27);

% gP1IS=econparams(28) no longer used****;

% dP1HH=econparams(29);zP1HH=econparams(30);kP1HH=econparams(31);

% mP1HH=econparams(32);nP1HH=econparams(33);

% dH1HH=econparams(34);zH1HH=econparams(35);kH1HH=econparams(36);

% mH1HH=econparams(37);nH1HH=econparams(38);

% dISHH=econparams(39);zISHH=econparams(40);kISHH=econparams(41);

% mISHH=econparams(42);nISHH=econparams(43);

% khat=econparams(44);theta=econparams(45);lambda=econparams(46);

% gRPP1=econparams(47);mP1=econparams(48);mH1=econparams(49);

% mHH=econparams(50);

% P1barmass=econparams(51);H1barmass=econparams(52);

% ISbarmass=econparams(53);dw=econparams(54);

% etaa=econparams(55);etab=econparams(56);etac=econparams(57);

% phi=econparams(58);

% idealpercapmass=econparams(59);

sol=[4.832248668616482e-001

1.357181037387526e-001

 1.000000000000000e+000 

60
 1.000000000000000e+000

 7.737088487049635e-001

 1.910770291618713e-001

 4.991249741897157e-002

 3.380538111199238e-002

 7.524098562805767e-001

 1.000000000000000e-003

 2.527165128459011e-001

 9.926234817539428e-001

 6.081040091094971e-001

 2.972103070158917e-001

 1.000000000000000e-003

 5.646817673172643e-001

 4.240782789383344e-001

 19.000000000000000e-001

 1.286507253686182e-001

 1.474467159083802e-002

 3.995283985198757e-001

 1.555235114499388e-001

 1.019919611353713e-001

 6.766772330024032e-001

 1.018295456166137e-003

 9.838862467656756e-003

 9.000000000000011e-001];

% sol(6)=1.25e-001;

%

% assign parameters

%

econparams(1:5)=sol(1:5);

econparams(6)=3*sol(6);econparams(7)=3*sol(7);

econparams(8:11)=sol(8:11);

econparams(12)=sol(6);

econparams(13)=sol(7);

econparams(14:17)=sol(12:15);

econparams(18)=sol(6);

econparams(19)=sol(7);

econparams(20)=sol(16);

dP1H1=sol(6);econparams(21:4:29)=[dP1H1/1000;dP1H1/1000;dP1H1/1000];eco

nparams(34:5:39)=[dP1H1/1000;dP1H1/1000];

eP1H1=sol(7);econparams(22:4:26)=[eP1H1;eP1H1];

fP1H1=sol(17);econparams(23:4:27)=[fP1H1/1;fP1H1/1];econparams(31:5:41)

=[fP1H1/1000;1/2*fP1H1/1000;1/2*fP1H1/1000];

econparams(24)=sol(18);econparams(28)=sol(19);

econparams(30:5:40)=[sol(20)/1000;sol(20)/1000;sol(20)/1000]*10000;

mH1HH=sol(21);econparams(32:5:42)=[1/2*mH1HH/1000;mH1HH/1000;1/2*mH1HH/

1000];

nISHH=sol(22);econparams(33:5:43)=[1/2*nISHH/1000;1/2*nISHH/1000;nISHH/

1000];

econparams(45)=sol(23);econparams(46)=sol(24);

econparams(48)=sol(25);econparams(49)=sol(26);econparams(50)=.22;

econparams(44)=0.3; %khat

econparams(47)=.09;%1.05e-2; %growth of P1 

61
econparams(51)=0;%1; %P1barmass

econparams(52)=.4;%0.125; %H1barmass

econparams(53)=0; %ISbarmass

econparams(54)=0.01*4.507354330528614e-004; %dw, second number in

% product is ic_econ(4)/numHH 

econparams(55)=9/6; %etaa

econparams(56)=5/6; %etab

econparams(57)=.1/6;%etac

econparams(58)=1e1; %phi 

% 

% initialize economic state [P1;H1;IS;HH]

% 

ic_econ=[1.550317582123858e+000

1.400236056583227e-001

 1.008970080362582e-001

 4.507354330528614e-001];

numHH=10; %initial number of humans 

econparams(59)=4.055586606327234e-003;%ic_econ(4)/numHH; %ideal per

% capita human mass 

% 

% calculate ecological equilibrium, without any of the domesticated or

% industrial agents

% 

ecolparams(8)=ecolparams(8)+.25*ecolparams(8);

ic_ecol=eq21apr2004V2(ecolparams); 

ecolparams(18)=0; %set RPIRP to zero, without changing original

% equilibrium. 

%comment this out to run ecological model itself 

% econparams=zeros(50,1); %uncomment these two lines to run the

% econparams(45:46,1)=[1;1]; %ecological model with no economics

% P1, P2, P3, H1, H2, H3, C1, C2, HH, IS, RP, IRP, P1H1massdeficit,

% P1ISmassdeficit,...

% P1HHmassdeficit, H1massdeficit, ISmassdeficit, numHH, percapmass

ic=[ic_econ(1);ic_ecol(1);ic_ecol(2);ic_econ(2);ic_ecol(3);ic_ecol(4);i

c_ecol(5);ic_ecol(6);ic_econ(4);...

ic_econ(3);ic_ecol(7);ic_ecol(8);0;0;0;0;0; numHH;

4.055586606327234e-003];%ic_econ(4)/numHH]; 

The ecological part of the model has its equilibrium calculated by the function 

eq21apr2004V2.m, given below. 

62
function out = eq21apr2004V2(params)

%

% eq21apr2004V2.m gives the equilibrium of the simple foodweb in

% equations3.nb, in C:\chrisp\EPA\projects\agent

% model\mnsimple_21_apr_2004

%

%

% First written 28 Jul 2004 CWP

%

% assign parameters

%

%

% params =

% [gRPP2;gP2H2;gP2H3;gRPP3;gP3H3;gH2C1;gH2C2;gH3C2;

% rIRPP2;rIRPP3;mP2;mP3;mH2;mH3;mC1;mC2;mIRPRP;gRPIRP;P2]; Here, P2 is

% a choice variable, controlling the total system mass

%

gRPP2=params(1);

gP2H2=params(2);gP2H3=params(3);gRPP3=params(4);gP3H3=params(5);

gH2C1=params(6);gH2C2=params(7);gH3C2=params(8);

rIRPP2=params(9);rIRPP3=params(10);

mP2=params(11);mP3=params(12);mH2=params(13);

mH3=params(14);mC1=params(15);mC2=params(16);

mIRPRP=params(17);RPIRP=params(18);P2=params(19); %gRPIRP not used here

%

% calculate equilibria as determined analytically in the Mathematica

% notebook

%

C1=(1/gH2C1)*(-mH2+gP2H2*P2-(1/gH3C2)*...

(gH2C2*(-mH3+gP2H3*P2+...

(gP3H3*((-gH2C2)*gP3H3*gRPP2*mC1*mIRPRP-...

gH3C2*gP2H2*gRPP3*mC1*mIRPRP+gH2C2*gP2H3*...

gRPP3*mC1*mIRPRP+gH2C1*gP3H3*gRPP2*mC2*...

mIRPRP-gH2C1*gP2H3*gRPP3*mC2*mIRPRP-...

gH2C1*gH3C2*gRPP3*mIRPRP*mP2+gH2C1*gH3C2*...

gRPP2*mIRPRP*mP3-gH2C2*gP3H3*gRPP2*mC1*P2*...

rIRPP2-gH3C2*gP2H2*gRPP3*mC1*P2*rIRPP2+...

gH2C2*gP2H3*gRPP3*mC1*P2*rIRPP2+...

gH2C1*gP3H3*gRPP2*mC2*P2*rIRPP2-...

gH2C1*gP2H3*gRPP3*mC2*P2*rIRPP2-...

gH2C1*gH3C2*gRPP3*mP2*P2*rIRPP2+...

gH2C1*gH3C2*gRPP2*mP3*P2*rIRPP2+...

gH2C1*gH3C2*gRPP3*rIRPP2*RPIRP-...

gH2C1*gH3C2*gRPP2*rIRPP3*RPIRP))/...

((gH2C2*gP3H3*gRPP2*mC1+gH3C2*gP2H2*gRPP3*mC1-...

gH2C2*gP2H3*gRPP3*mC1-gH2C1*gP3H3*gRPP2*mC2+...

gH2C1*gP2H3*gRPP3*mC2+gH2C1*gH3C2*gRPP3*mP2-...

gH2C1*gH3C2*gRPP2*mP3)*rIRPP3))));

C2=(1/gH3C2)*(-mH3+gP2H3*P2+...

(gP3H3*((-gH2C2)*gP3H3*gRPP2*mC1*mIRPRP-...

gH3C2*gP2H2*gRPP3*mC1*mIRPRP+gH2C2*gP2H3*gRPP3*...

63
mC1*mIRPRP+gH2C1*gP3H3*gRPP2*mC2*mIRPRP-...

gH2C1*gP2H3*gRPP3*mC2*mIRPRP-gH2C1*gH3C2*gRPP3*...

mIRPRP*mP2+gH2C1*gH3C2*gRPP2*mIRPRP*mP3-...

gH2C2*gP3H3*gRPP2*mC1*P2*rIRPP2-...

gH3C2*gP2H2*gRPP3*mC1*P2*rIRPP2+...

gH2C2*gP2H3*gRPP3*mC1*P2*rIRPP2+...

gH2C1*gP3H3*gRPP2*mC2*P2*rIRPP2-...

gH2C1*gP2H3*gRPP3*mC2*P2*rIRPP2-...

gH2C1*gH3C2*gRPP3*mP2*P2*rIRPP2+...

gH2C1*gH3C2*gRPP2*mP3*P2*rIRPP2+...

gH2C1*gH3C2*gRPP3*rIRPP2*RPIRP-...

gH2C1*gH3C2*gRPP2*rIRPP3*RPIRP))/...

((gH2C2*gP3H3*gRPP2*mC1+gH3C2*gP2H2*gRPP3*mC1-...

gH2C2*gP2H3*gRPP3*mC1-gH2C1*gP3H3*gRPP2*mC2+...

gH2C1*gP2H3*gRPP3*mC2+gH2C1*gH3C2*gRPP3*mP2-...

gH2C1*gH3C2*gRPP2*mP3)*rIRPP3)); 

H2=mC1/gH2C1; 

H3=((-gH2C2)*mC1 + gH2C1*mC2)/(gH2C1*gH3C2); 

P3=((-gH2C2)*gP3H3*gRPP2*mC1*mIRPRP-gH3C2*gP2H2*gRPP3*mC1*...

mIRPRP+gH2C2*gP2H3*gRPP3*mC1*mIRPRP+...

gH2C1*gP3H3*gRPP2*mC2*mIRPRP-gH2C1*gP2H3*gRPP3*mC2*...

mIRPRP-gH2C1*gH3C2*gRPP3*mIRPRP*mP2+...

gH2C1*gH3C2*gRPP2*mIRPRP*mP3-gH2C2*gP3H3*gRPP2*mC1*...

P2*rIRPP2-gH3C2*gP2H2*gRPP3*mC1*P2*rIRPP2+...

gH2C2*gP2H3*gRPP3*mC1*P2*rIRPP2+...

gH2C1*gP3H3*gRPP2*mC2*P2*rIRPP2-...

gH2C1*gP2H3*gRPP3*mC2*P2*rIRPP2-...

gH2C1*gH3C2*gRPP3*mP2*P2*rIRPP2+...

gH2C1*gH3C2*gRPP2*mP3*P2*rIRPP2+...

gH2C1*gH3C2*gRPP3*rIRPP2*RPIRP-...

gH2C1*gH3C2*gRPP2*rIRPP3*RPIRP)/...

((gH2C2*gP3H3*gRPP2*mC1+gH3C2*gP2H2*gRPP3*mC1-...

gH2C2*gP2H3*gRPP3*mC1-gH2C1*gP3H3*gRPP2*mC2+...

gH2C1*gP2H3*gRPP3*mC2+gH2C1*gH3C2*gRPP3*mP2-...

gH2C1*gH3C2*gRPP2*mP3)*rIRPP3); 

RP=((-gH2C2)*gP3H3*mC1*rIRPP2+gH2C1*gP3H3*mC2*rIRPP2+...

gH2C1*gH3C2*mP3*rIRPP2-gH3C2*gP2H2*mC1*rIRPP3+...

gH2C2*gP2H3*mC1*rIRPP3-gH2C1*gP2H3*mC2*rIRPP3-...

gH2C1*gH3C2*mP2*rIRPP3)/(gH2C1*gH3C2*...

(gRPP3*rIRPP2-gRPP2*rIRPP3)); 

IRP=((-gH2C2)*gP3H3*gRPP2*mC1-gH3C2*gP2H2*gRPP3*mC1+...

gH2C2*gP2H3*gRPP3*mC1+gH2C1*gP3H3*gRPP2*mC2-...

gH2C1*gP2H3*gRPP3*mC2-gH2C1*gH3C2*gRPP3*mP2+...

gH2C1*gH3C2*gRPP2*mP3)/(gH2C1*gH3C2*...

((-gRPP3)*rIRPP2+gRPP2*rIRPP3)); 

out=[P2;P3;H2;H3;C1;C2;RP;IRP]; 

S- Function, core of model (file: price_settingSv8.m) 

64
function [sys,y0,str,ts] =

price_settingSv8(t,y,u,flag,ecolparams,econparams,...

iconditions,MH_updown,Lamb_updown,Thet_updown, tlo, thi)

% price_settingSv8.m is based on price_settingSv7.m. corrects v7 

% See summary summtoaustriav5.doc

% 

% 

% 

% 

% 

% 18 Jan 2005 CWP 

% 

%

% The following outlines the general structure of an S-function.

%

switch flag,

 %%%%%%%%%%%%%%%%%%

 % Initialization %

 %%%%%%%%%%%%%%%%%%

 case 0,

[sys,y0,str,ts]=mdlInitializeSizes(iconditions,ecolparams,econparams);

 %%%%%%%%%%%%%%%

 % Derivatives %

 %%%%%%%%%%%%%%%

 case 1,

sys=mdlDerivatives(t,y,u);

 %%%%%%%%%%

 % Update %

%%%%%%%%%%

 case 2,

sys=mdlUpdate(t,y,u,ecolparams,econparams,MH_updown,Lamb_updown,Thet_up

down, tlo, thi);

 %%%%%%%%%%%

 % Outputs %

%%%%%%%%%%%

 case 3,

sys=mdlOutputs(t,y,u);

 %%%%%%%%%%%%%%%%%%%%%%%

 % GetTimeOfNextVarHit %

 %%%%%%%%%%%%%%%%%%%%%%%

 case 4,

sys=mdlGetTimeOfNextVarHit(t,y,u);

 %%%%%%%%%%%%%

 % Terminate %

 %%%%%%%%%%%%%

65
 case 9,

sys=mdlTerminate(t,y,u);

 %%%%%%%%%%%%%%%%%%%%

 % Unexpected flags %

%%%%%%%%%%%%%%%%%%%%

 otherwise

 error(['Unhandled flag = ',num2str(flag)]); 

end

% end sfuntmpl

%

%======================================================================

% mdlInitializeSizes

% Return the sizes, initial conditions, and sample times for the S-

% function.

%======================================================================

%

function 

[sys,y0,str,ts]=mdlInitializeSizes(iconditions,ecolparams,econparams)

global RPP1 P1H2 P1IS P1H1 P1HH RPIS ISIRP P2H1 H1C1 H1HH pP1 pH1 pIS W

ISHHflow

%

% call simsizes for a sizes structure, fill it in and convert it to a

% sizes array.

%

sizes = simsizes;

sizes.NumContStates = 0;

sizes.NumDiscStates = 19;

sizes.NumOutputs = 33;

sizes.NumInputs = 0;

sizes.DirFeedthrough = 0;

sizes.NumSampleTimes = 1; % at least one sample time is needed

sys = simsizes(sizes);

%

% assign initial conditions

%

y0=iconditions;

%

% str is always an empty matrix

%

str = [];

%

% initialize the array of sample times

%

ts = [1 0];

66
%

% Do the following to find initial transfer flows for the global

% variables

%

belownoreproduction=1e-4; %level below which the natural ecosystem

% elements do not reproduce

% assign parameters

%

%

% ecolparams =

% [gRPP2;gP2H2;gP2H3;gRPP3;gP3H3;gH2C1;gH2C2;gH3C2;

% rIRPP2;rIRPP3;mP2;mP3;mH2;mH3;mC1;mC2;mIRPRP;RPIRP;P2]; Here, P2 is a

% choice variable, controlling the total system mass

%

gRPP2=ecolparams(1);

gP2H2=ecolparams(2);gP2H3=ecolparams(3);gRPP3=ecolparams(4);gP3H3=ecolp

arams(5);

gH2C1=ecolparams(6);gH2C2=ecolparams(7);gH3C2=ecolparams(8);

rIRPP2=ecolparams(9);rIRPP3=ecolparams(10);

mP2=ecolparams(11);mP3=ecolparams(12);mH2=ecolparams(13);

mH3=ecolparams(14);mC1=ecolparams(15);mC2=ecolparams(16);

mIRPRP=ecolparams(17);RPIRP=ecolparams(18); %this last constant changed

% 2 Aug 2004

gP1H2=ecolparams(20);gH1C1=ecolparams(21); %13 Aug 2004 natural

% predation of P1 and H1

%

% econparams are economic parameters:

%

aw=econparams(1);cw=econparams(2);

aP1=econparams(3);bP1=econparams(4);cP1=econparams(5);

aP1p=econparams(6);bP1p=econparams(7);cP1p=econparams(8);

aH1=econparams(9);bH1=econparams(10);cH1=econparams(11);

aH1p=econparams(12);bH1p=econparams(13);cH1p=econparams(14);

aIS=econparams(15);bIS=econparams(16);cIS=econparams(17);

aISp=econparams(18);bISp=econparams(19);cISp=econparams(20);

dP1H1=econparams(21);eP1H1=econparams(22);fP1H1=econparams(23);gP1H1=ec

onparams(24);

% dP1IS=econparams(25);eP1IS=econparams(26);fP1IS=econparams(27);

% gP1IS=econparams(28);

dP1HH=econparams(29);zP1HH=econparams(30);kP1HH=econparams(31);mP1HH=ec

onparams(32);nP1HH=econparams(33);

dH1HH=econparams(34);zH1HH=econparams(35);kH1HH=econparams(36);mH1HH=ec

onparams(37);nH1HH=econparams(38);

dISHH=econparams(39);zISHH=econparams(40);kISHH=econparams(41);mISHH=ec

onparams(42);nISHH=econparams(43);

khat=econparams(44);theta=econparams(45);lambda=econparams(46);

gRPP1=econparams(47);mP1=econparams(48);mH1=econparams(49);mHH=econpara

ms(50);

P1bar=econparams(51);H1bar=econparams(52);ISbar=econparams(53);dw=econp

arams(54);

etaa=econparams(55);etab=econparams(56);etac=econparams(57);phi=econpar

ams(58);

67
idealpercapmass=econparams(59);

%

% assign state

%

P1=y0(1);P2=y0(2);P3=y0(3);

H1=y0(4);H2=y0(5);H3=y0(6);

C1=y0(7);C2=y0(8);

HH=y0(9);ISmass=y0(10); %ISmassdeficit keeps track of deficits

RP=y0(11);

IRP=y0(12);

P1H1massdeficit=y0(13);P1ISmassdeficit=y0(14);P1HHmassdeficit=y0(15);

H1massdeficit=y0(16);ISmassdeficit=y0(17);

numHH=y0(18);percapmass=y0(19);

P1massdeficit=P1H1massdeficit+P1ISmassdeficit+P1HHmassdeficit;

rIRPP2=rIRPP2*(10^2/(10^2+IRP^2));

rIRPP3=rIRPP3*(10^2/(10^2+IRP^2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% % 

% Economics % 

% % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% Industrial sector sets the wage rate

% 

W=max(aw+cw*(ISbar-(ISmassdeficit+ISmass))/(theta+lambda)-dw*numHH,0); 

% 

% Based on the Wage, industries set prices and their production (how

% much they would like to produce to maximize their profits based on 

% their assumption as to what the demand for their products will be).

% Here a linear functional form for the supply is assumed.

% 

if P1==0

 pP1=0;

P1production=0;

else

 pP1=max(aP1+bP1*W-cP1*((P1massdeficit+P1)-P1bar),0);

P1production=max(aP1p-bP1p*W-cP1p*((P1massdeficit+P1)-P1bar),0);

end 

if H1==0

 pH1=0;

H1production=0;

else

 pH1=max(aH1+bH1*W-cH1*((H1massdeficit+H1)-H1bar),0);

H1production=max(aH1p-bH1p*W-cH1p*((H1massdeficit+H1)-H1bar),0);

end 

pIS=max(aIS+bIS*W+cIS*(ISbar-(ISmassdeficit+ISmass))/(theta+lambda),0);

ISproduction=max(aISp-bISp*W+cISp*(ISbar-

(ISmassdeficit+ISmass))/(theta+lambda),0); 

68
if HH==0|numHH<2

pIS=0;

ISproduction=0;

end 

% 

% Next, a calculation is made of how much each industry is going to

% demand of its suppliers. The following are in units of mass, unless

% otherwise noted 

% 

if H1==0|HH==0|numHH<2

P1H1demand=0;

P2H1=0;

else

 P1H1demand=max(dP1H1-eP1H1*W-fP1H1*pP1-gP1H1*((H1massdeficit+H1)-

H1bar),0);

P2H1=khat;

end 

%

% expressions for P1HH, H1HH and ISHH reflect constraint on human

% spending.

% These expressions were determined in MATHEMATICA under the file

% C:\chrisp\EPA\projects\agent model\mnsimple_21_apr_2004

% solve_humans.nb

% These are per capita, so must be multiplied by population later

%

P1HHdemand=max((1/(-1+zP1HH+zH1HH+zISHH))*(-dP1HH-mP1HH*pH1-

nP1HH*pIS+...

kP1HH*pP1-dH1HH*zP1HH-dISHH*zP1HH+mH1HH*pH1*zP1HH-...

mISHH*pH1*zP1HH-nH1HH*pIS*zP1HH+nISHH*pIS*zP1HH-...

kH1HH*pP1*zP1HH-kISHH*pP1*zP1HH+dP1HH*zH1HH+...

mP1HH*pH1*zH1HH+nP1HH*pIS*zH1HH-kP1HH*pP1*zH1HH+...

dP1HH*zISHH+mP1HH*pH1*zISHH+nP1HH*pIS*zISHH-kP1HH*pP1*zISHH),0);

H1HHdemand=max((1/(-1+zP1HH+zH1HH+zISHH))*(-dH1HH+mH1HH*pH1-nH1HH*pIS-

...

kH1HH*pP1+dH1HH*zP1HH-mH1HH*pH1*zP1HH+nH1HH*pIS*zP1HH+...

kH1HH*pP1*zP1HH-dISHH*zH1HH-dP1HH*zH1HH-mISHH*pH1*zH1HH-...

mP1HH*pH1*zH1HH+nISHH*pIS*zH1HH-nP1HH*pIS*zH1HH-...

kISHH*pP1*zH1HH+kP1HH*pP1*zH1HH+dH1HH*zISHH-...

mH1HH*pH1*zISHH+nH1HH*pIS*zISHH+kH1HH*pP1*zISHH),0);

ISHHdemand=max(-((dISHH+mISHH*pH1-nISHH*pIS+kISHH*pP1-dISHH*zP1HHmISHH*pH1*zP1HH+...

nISHH*pIS*zP1HH-kISHH*pP1*zP1HH-dISHH*zH1HHmISHH*pH1*zH1HH+nISHH*pIS*zH1HH-...

kISHH*pP1*zH1HH+dH1HH*zISHH+dP1HH*zISHHmH1HH*pH1*zISHH+mP1HH*pH1*zISHH+...

nH1HH*pIS*zISHH+nP1HH*pIS*zISHH+kH1HH*pP1*zISHH-kP1HH*pP1*zISHH)/...

(-1+zP1HH+zH1HH+zISHH)),0); %in units of units

% corrected 9 Aug 2004

if HH==0|numHH<2

P1HHdemand=0; 

69
 H1HHdemand=0;

ISHHdemand=0;

end 

% 

% The flows that involve labor to keep the wild

% from taking domestics, namely P1H2 and H1C2 must then be calculated.

% 

if P1==0|H2==0

P1H2=0;

else

 P1H2=max((gRPP1*P1*RP-mP1*P1-P1production),0);

end 

if HH==0|numHH<2

P1H2=gP1H2*P1*H2;

end 

if H1==0|C1==0

H1C1=0;

else

 H1C1=max((P1H1demand+P2H1-mH1*H1-H1production),0);

end 

if HH==0|numHH<2

H1C1=gH1C1*H1*C1;

end 

P1ISdemand=theta*ISproduction;

RPISdemand=lambda*ISproduction; 

% 

% calculate all but next state, according to system

% equations (pricesettingequations3.doc, and Whitmore's paper

% C:\chrisp\EPA\projects\agent model\mnsimple_21_apr_2004\whitdocs\12

% cell imp comp EPA 5-04-04.doc)

% 

% 

% P1 

% 

P1RP=max(mP1*P1,0);RPP1=max(gRPP1*P1*RP,0);

P1H1=P1H1demand;

P1IS=P1ISdemand;

P1HH=P1HHdemand*numHH;

if P1+RPP1-P1RP-P1H2-P1H1-P1HH-P1IS<0 %if statement to deal with going

% negative

if P1+RPP1-P1RP<0

 P1RP=P1+RPP1;

P1H2=0;P1H1=0;P1HH=0;P1IS=0;

else

 totP1demand=P1H2+P1H1+P1HH+P1IS;

P1avail=P1+RPP1-P1RP;

P1H2=P1avail*P1H2/totP1demand;

70
 P1H1=P1avail*P1H1/totP1demand;

P1HH=P1avail*P1HH/totP1demand;

P1IS=P1avail-(P1H2+P1H1+P1HH);

end 

else

 if P1massdeficit<0 %if there is an accumulated deficit between 

% demand for P1 try to

%make this up if there is extra stock

P1surplus=min(P1+RPP1-P1RP-P1H2-P1H1-P1HH-P1IS,-P1massdeficit);

%only what you need to make up deficit

P1H1=P1H1+P1surplus*P1H1massdeficit/P1massdeficit;

P1IS=P1IS+P1surplus*P1ISmassdeficit/P1massdeficit;

P1HH=P1HH+P1surplus*P1HHmassdeficit/P1massdeficit;

end 

end

P2H2=gP2H2*P2*H2;P2H3=gP2H3*P2*H3;P2RP=max(mP2*P2,0);RPP2=max(gRPP2*RP*

P2,0);

IRPP2=max(rIRPP2*P2*IRP,0);

P3RP=max(mP3*P3,0);P3H3=gP3H3*P3*H3;RPP3=max(gRPP3*RP*P3,0);

IRPP3=max(rIRPP3*P3*IRP,0);

if IRP<=0

 IRPP2=0;

IRPP3=0;

elseif IRP-IRPP2-IRPP3-max(IRP*mIRPRP,0)+RPIRP<0

if P2~=0

 IRPP2=rIRPP2*(IRP-max(IRP*mIRPRP,0)+RPIRP)/(rIRPP2+rIRPP3);

end

 if P3~=0

 IRPP3=rIRPP3*(IRP-max(IRP*mIRPRP,0)+RPIRP)/(rIRPP2+rIRPP3);

end 

end 

if P2+IRPP2+RPP2-P2RP-P2H2-P2H3-P2H1<belownoreproduction

if P2+IRPP2+RPP2-P2RP<belownoreproduction

P2RP=P2+IRPP2+RPP2;

P2H2=0;P2H3=0;P2H1=0;

else

 totP2demand=P2H2+P2H3+P2H1;

P2avail=P2+IRPP2+RPP2-P2RP;

P2H2=P2H2*P2avail/totP2demand;

P2H3=P2H3*P2avail/totP2demand;

P2H1=P2avail-(P2H2+P2H3);

end

end 

if P3+IRPP3+RPP3-P3RP-P3H3<belownoreproduction

if P3+IRPP3+RPP3-P3RP<belownoreproduction

P3RP=P3+IRPP3+RPP3;

P3H3=0;

else

 totP3demand=P3H3;

P3avail=P3+IRPP3+RPP3-P3RP;

P3H3=P3H3*P3avail/totP3demand;

end 

end 

71
H1RP=max(mH1*H1,0);

H1HH=H1HHdemand*numHH;

if H1+P1H1+P2H1-H1RP-H1C1-H1HH<0

 if H1+P1H1+P2H1-H1RP<0

 H1RP=H1+P1H1+P2H1;

H1C1=0;H1HH=0;

else

 totH1demand=H1C1+H1HH;

H1avail=H1+P1H1+P2H1-H1RP;

H1C1=H1avail*H1C1/totH1demand;

H1HH=H1avail-H1C1;

end 

else

 if H1massdeficit<0

 H1HH=H1HH+min(H1+P1H1+P2H1-H1RP-H1C1-H1HH,-H1massdeficit);

%only what you need to make up deficit

end 

end 

H2C1=gH2C1*C1*H2;H2C2=gH2C2*H2*C2;H2RP=max(mH2*H2,0);

if H2+P1H2+P2H2-H2RP-H2C1-H2C2<belownoreproduction

if H2+P1H2+P2H2-H2RP<belownoreproduction

H2RP=H2+P1H2+P2H2;

H2C1=0;H2C2=0;

else

 totH2demand=H2C1+H2C2;

H2avail=H2+P1H2+P2H2-H2RP;

H2C1=H2C1*H2avail/totH2demand;

H2C2=H2avail-H2C1;

end 

end 

H3RP=max(mH3*H3,0);H3C2=gH3C2*H3*C2;

if H3+P2H3+P3H3-H3RP-H3C2<belownoreproduction

if H3+P2H3+P3H3-H3RP<belownoreproduction

H3RP=H3+P2H3+P3H3;

H3C2=0;

else

 totH3demand=H3C2;

H3avail=H3+P2H3+P3H3-H3RP;

H3C2=H3C2*H3avail/totH3demand;

end 

end 

C1RP=max(mC1*C1,0);

if C1+H1C1+H2C1-C1RP<belownoreproduction

C1RP=C1+H1C1+H2C1;

end 

C2RP=max(mC2*C2,0);

if C2+H2C2+H3C2-C2RP<belownoreproduction

C2RP=C2+H2C2+H3C2;

end 

HHRP=mHH*numHH*percapmass; 

72
IRPRP=max(IRP*mIRPRP,0);

RPIS=min(lambda*P1IS/theta,RPISdemand);

stockRP=RP+P1RP+P2RP+P3RP+H1RP+H2RP+H3RP+C1RP+C2RP+HHRP+IRPRP;

if stockRP<0

 stockRP=0;

end 

if stockRP-(RPP1+RPP2+RPP3)-RPIRP-RPIS<=0&RPIRP==0 %this line changed

% 26 Aug 2004

RPdemand=RPP1+RPP2+RPP3+RPIS;

RPP1=RPP1*stockRP/RPdemand;

RPP2=RPP2*stockRP/RPdemand;

RPP3=RPP3*stockRP/RPdemand;

if RPIS~=0

 RPIS=stockRP-(RPP1+RPP2+RPP3);

else

 RPIS=0;

end 

end 

P1IS=min(theta*RPIS/lambda,P1IS);

% t 

% stockRP 

% nextRP = stockRP-(RPP1+RPP2+RPP3)-RPIRP-RPIS 

ISHHflow=max(0,(theta+lambda)*ISHHdemand*numHH);

ISIRP=ISHHflow;

if ISmass+P1IS+RPIS-ISIRP<=0

 ISIRP=ISmass+P1IS+RPIS;

else

 if ISmassdeficit<0&numHH>=2 %if there is an accumulated deficit 

% between demand for IS by the HH and IS

%supplied, try to make this up if there is extra stock

ISIRP=ISIRP+min(ISmass+P1IS+RPIS-ISIRP,-ISmassdeficit); %only

% what you need to make up deficit

end 

end 

if (P1HHdemand+H1HHdemand+ISHHflow)==0

weightedprice=0;

percapbirths=0;

elseif (pP1*P1HHdemand+pH1*H1HHdemand+pIS*ISHHflow)==0

weightedprice=0;

percapbirths=0;

else 

%weightedprice=(pP1*P1HHdemand+pH1*H1HHdemand+pIS*ISHHflow)/(P1HHdemand

% +H1HHdemand+ISHHflow);

weightedprice=(pP1*P1HH+pH1*H1HH+pIS*ISIRP)/(P1HH+H1HH+ISIRP);

%percapbirths=max(etaa-

% etab*W/weightedprice+etac*(W/weightedprice)^2,0);

percapbirths=max(etaa-etab*sqrt(W/weightedprice),0);

% if etab<=2*etac*(W/weightedprice)|etaa<=etac*(W/weightedprice)^2

% error('pricesettingmodel:growtherror','Percapbirths function

% parameters give bad function form')

% end 

73
end 

% end mdlInitializeSizes 

%

%======================================================================

% mdlDerivatives

% Return the derivatives for the continuous states.

%======================================================================

%

function sys=mdlDerivatives(t,y,u)

sys=[];

% end mdlDerivatives

%

%======================================================================

% mdlUpdate

% Handle discrete state updates, sample time hits, and major time step

% requirements.

%======================================================================

%

function 

sys=mdlUpdate(t,y,u,ecolparams,econparams,MH_updown,Lamb_updown,Thet_up

down, tlo, thi)

global RPP1 P1H2 P1IS P1H1 P1HH RPIS ISIRP P2H1 H1C1 H1HH pP1 pH1 pIS W

ISHHflow

belownoreproduction=1e-4; %level below which the natural ecosystem

% elements do not reproduce

% assign parameters

%

%

% ecolparams =

% [gRPP2;gP2H2;gP2H3;gRPP3;gP3H3;gH2C1;gH2C2;gH3C2;

% rIRPP2;rIRPP3;mP2;mP3;mH2;mH3;mC1;mC2;mIRPRP;RPIRP;P2]; Here, P2 is a

% choice variable, controlling the total system mass %gRPIRP changed 2

% Aug 2004 to RPIRP

%

gRPP2=ecolparams(1);

gP2H2=ecolparams(2);gP2H3=ecolparams(3);gRPP3=ecolparams(4);gP3H3=ecolp

arams(5);

gH2C1=ecolparams(6);gH2C2=ecolparams(7);gH3C2=ecolparams(8);

rIRPP2=ecolparams(9);rIRPP3=ecolparams(10);

mP2=ecolparams(11);mP3=ecolparams(12);mH2=ecolparams(13);

74
mH3=ecolparams(14);mC1=ecolparams(15);mC2=ecolparams(16);

mIRPRP=ecolparams(17);RPIRP=ecolparams(18);

gP1H2=ecolparams(20);gH1C1=ecolparams(21); %13 Aug 2004 natural

predation of P1 and H1

%

% econparams are economic parameters:

%

aw=econparams(1);cw=econparams(2);

aP1=econparams(3);bP1=econparams(4);cP1=econparams(5);

aP1p=econparams(6);bP1p=econparams(7);cP1p=econparams(8);

aH1=econparams(9);bH1=econparams(10);cH1=econparams(11);

aH1p=econparams(12);bH1p=econparams(13);cH1p=econparams(14);

aIS=econparams(15);bIS=econparams(16);cIS=econparams(17);

aISp=econparams(18);bISp=econparams(19);cISp=econparams(20);

dP1H1=econparams(21);eP1H1=econparams(22);fP1H1=econparams(23);gP1H1=ec

onparams(24);

%dP1IS=econparams(25);eP1IS=econparams(26);fP1IS=econparams(27);gP1IS=e

% conparams(28);

dP1HH=econparams(29);zP1HH=econparams(30);kP1HH=econparams(31);mP1HH=ec

onparams(32);nP1HH=econparams(33);

dH1HH=econparams(34);zH1HH=econparams(35);kH1HH=econparams(36);mH1HH=ec

onparams(37);nH1HH=econparams(38);

dISHH=econparams(39);zISHH=econparams(40);kISHH=econparams(41);mISHH=ec

onparams(42);nISHH=econparams(43);

khat=econparams(44);theta=econparams(45);lambda=econparams(46);

gRPP1=econparams(47);mP1=econparams(48);mH1=econparams(49);mHH=econpara

ms(50);

P1bar=econparams(51);H1bar=econparams(52);ISbar=econparams(53);dw=econp

arams(54);

etaa=econparams(55);etab=econparams(56);etac=econparams(57);phi=econpar

ams(58);

idealpercapmass=econparams(59);

% gRPP1=gRPP1*(1+1/100*sin(2*pi*t/12));

% gRPP2=gRPP2*(1+1/100*sin(2*pi*t/12));

% gRPP3=gRPP3*(1+1/100*sin(2*pi*t/12));

lambdaz=lambda;thetaz=theta;mHHz=mHH;

if Lamb_updown~=1

if Lamb_updown==2

if t>=tlo&t<=thi

 lambda=lambdaz+.001*(t-tlo);

elseif t>thi

 lambda=lambdaz+.001*(thi-tlo);

end

 else

 if t>=tlo

 lambda=max(lambdaz-.001*(t-tlo),lambda/10);

end

 end 

end 

75
if Thet_updown~=1

if Thet_updown==2

if t>=tlo

 theta=thetaz+.001*(t-tlo);

end

 else

 if t>=tlo

 theta=max(thetaz-.001*(t-tlo),theta/10);

end

 end 

end 

if MH_updown~=1

if MH_updown==2

if t>=tlo&t<=thi

 mHH=min(mHH*2,mHHz+.001*(t-tlo));

elseif t>thi

 mHH=min(mHH*2,mHHz+.001*(thi-tlo));

end

 else

 if t>=tlo

 mHH=max(mHH/10,mHHz-.001*(t-tlo));

end

 end 

end 

%

% assign state

%

P1=y(1);P2=y(2);P3=y(3);

H1=y(4);H2=y(5);H3=y(6);

C1=y(7);C2=y(8);

HH=y(9);ISmass=y(10);

RP=y(11);

IRP=y(12);

P1H1massdeficit=y(13);P1ISmassdeficit=y(14);P1HHmassdeficit=y(15);

H1massdeficit=y(16);ISmassdeficit=y(17);

numHH=y(18);percapmass=y(19);

P1massdeficit=P1H1massdeficit+P1ISmassdeficit+P1HHmassdeficit;

rIRPP2=rIRPP2*(10^2/(10^2+IRP^2));

rIRPP3=rIRPP3*(10^2/(10^2+IRP^2)); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% % 

% Economics % 

% % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

76
% Industrial sector sets the wage rate

% 

W=max(aw+cw*(ISbar-(ISmassdeficit+ISmass))/(theta+lambda)-dw*numHH,0); 

% 

% Based on the Wage, industries set prices and their production (how

% much they would like to produce to maximize their profits based on 

% their assumption as to what the demand for their products will be).

% Here, a linear functional form is assumed for the supply. if P1==0

 pP1=0;

P1production=0;

else

 pP1=max(aP1+bP1*W-cP1*((P1massdeficit+P1)-P1bar),0);

P1production=max(aP1p-bP1p*W-cP1p*((P1massdeficit+P1)-P1bar),0);

end 

if H1==0

 pH1=0;

H1production=0;

else

 pH1=max(aH1+bH1*W-cH1*((H1massdeficit+H1)-H1bar),0);

H1production=max(aH1p-bH1p*W-cH1p*((H1massdeficit+H1)-H1bar),0);

end

pIS=max(aIS+bIS*W+cIS*(ISbar-(ISmassdeficit+ISmass))/(theta+lambda),0);

ISproduction=max(aISp-bISp*W+cISp*(ISbar-

(ISmassdeficit+ISmass))/(theta+lambda),0);

if HH==0|numHH<2

pIS=0;

ISproduction=0;

end 

% 

% Next, how much each industry is going to demand of its

% suppliers is calculated. The following are in units of mass, unless

% otherwise noted 

% 

if H1==0|HH==0|numHH<2

P1H1demand=0;

P2H1=0;

else

 P1H1demand=max(dP1H1-eP1H1*W-fP1H1*pP1-gP1H1*((H1massdeficit+H1)-

H1bar),0);

P2H1=khat;

end 

% 

% expressions for P1HH, H1HH and ISHH reflect constraint on human

% spending.

% These expressions were determined in MATHEMATICA under the file

% C:\chrisp\EPA\projects\price_setting\mnsimple_21_apr_2004

% solve_humans.nb 

% These are per capita, so must be multiplied by population later 

77
%

P1HHdemand=max((1/(-1+zP1HH+zH1HH+zISHH))*(-dP1HH-mP1HH*pH1-

nP1HH*pIS+...

kP1HH*pP1-dH1HH*zP1HH-dISHH*zP1HH+mH1HH*pH1*zP1HH-...

mISHH*pH1*zP1HH-nH1HH*pIS*zP1HH+nISHH*pIS*zP1HH-...

kH1HH*pP1*zP1HH-kISHH*pP1*zP1HH+dP1HH*zH1HH+...

mP1HH*pH1*zH1HH+nP1HH*pIS*zH1HH-kP1HH*pP1*zH1HH+...

dP1HH*zISHH+mP1HH*pH1*zISHH+nP1HH*pIS*zISHH-kP1HH*pP1*zISHH),0);

H1HHdemand=max((1/(-1+zP1HH+zH1HH+zISHH))*(-dH1HH+mH1HH*pH1-nH1HH*pIS-

...

kH1HH*pP1+dH1HH*zP1HH-mH1HH*pH1*zP1HH+nH1HH*pIS*zP1HH+...

kH1HH*pP1*zP1HH-dISHH*zH1HH-dP1HH*zH1HH-mISHH*pH1*zH1HH-...

mP1HH*pH1*zH1HH+nISHH*pIS*zH1HH-nP1HH*pIS*zH1HH-...

kISHH*pP1*zH1HH+kP1HH*pP1*zH1HH+dH1HH*zISHH-...

mH1HH*pH1*zISHH+nH1HH*pIS*zISHH+kH1HH*pP1*zISHH),0);

ISHHdemand=max(-((dISHH+mISHH*pH1-nISHH*pIS+kISHH*pP1-dISHH*zP1HHmISHH*pH1*zP1HH+...

nISHH*pIS*zP1HH-kISHH*pP1*zP1HH-dISHH*zH1HHmISHH*pH1*zH1HH+nISHH*pIS*zH1HH-...

kISHH*pP1*zH1HH+dH1HH*zISHH+dP1HH*zISHHmH1HH*pH1*zISHH+mP1HH*pH1*zISHH+...

nH1HH*pIS*zISHH+nP1HH*pIS*zISHH+kH1HH*pP1*zISHH-kP1HH*pP1*zISHH)/...

(-1+zP1HH+zH1HH+zISHH)),0); %in units of units

% corrected 9 Aug 2004 

if HH==0|numHH<2

ISHHdemand=0;

P1HHdemand=0;

H1HHdemand=0;

end 

% ISHHdemand+P1HHdemand+H1HHdemand 

% P1HHdemand 

% 

% The flows that involve labor to keep the wild

% from taking domestics, namely P1H2 and H1C2 must then be calculated.

% 

if P1==0|H2==0

P1H2=0;

else

 P1H2=max((gRPP1*P1*RP-mP1*P1-P1production),0);

end 

if HH==0|numHH<2

P1H2=gP1H2*P1*H2;

end 

if H1==0|C1==0

H1C1=0;

else

 H1C1=max((P1H1demand+P2H1-mH1*H1-H1production),0);

end 

78
if HH==0|numHH<2

H1C1=gH1C1*H1*C1;

end 

% 

% the ISproduction is checked again below as well, after checks for

% realistic mass transfers 

% 

if HH==0|numHH<2

ISproduction=0.0;

end 

P1ISdemand=theta*ISproduction;

RPISdemand=lambda*ISproduction; 

% 

% calculate next state, according to system

% equations (pricesettingequations3.doc, and Whitmore's paper

% C:\chrisp\EPA\projects\agent model\mnsimple_21_apr_2004\whitdocs\12

% cell imp comp EPA 5-04-04.doc)

% Here, check to see that these transfers won't violate conservation of

% mass 

% 

% P1 

% 

P1RP=max(mP1*P1,0);RPP1=max(gRPP1*P1*RP,0);

P1H1=P1H1demand;

P1IS=P1ISdemand;

P1HH=P1HHdemand*numHH;

if P1+RPP1-P1RP-P1H2-P1H1-P1HH-P1IS<0 %if statement to deal with going

% negative

if P1+RPP1-P1RP<0

 P1RP=P1+RPP1;

P1H2=0;P1H1=0;P1HH=0;P1IS=0;

else

 totP1demand=P1H2+P1H1+P1HH+P1IS;

P1avail=P1+RPP1-P1RP;

P1H2=P1avail*P1H2/totP1demand;

P1H1=P1avail*P1H1/totP1demand;

P1HH=P1avail*P1HH/totP1demand;

P1IS=P1avail-(P1H2+P1H1+P1HH);

end 

else

 if P1massdeficit<0 %if there is an accumulated deficit between 

% demand for P1 try to

%make this up if there is extra stock

P1surplus=min(P1+RPP1-P1RP-P1H2-P1H1-P1HH-P1IS,-P1massdeficit);

%only what you need to make up deficit

P1H1=P1H1+P1surplus*P1H1massdeficit/P1massdeficit;

P1IS=P1IS+P1surplus*P1ISmassdeficit/P1massdeficit;

P1HH=P1HH+P1surplus*P1HHmassdeficit/P1massdeficit;

end 

end 

% 

% P2 

79
% 

P2H2=gP2H2*P2*H2;P2H3=gP2H3*P2*H3;P2RP=max(mP2*P2,0);RPP2=max(gRPP2*RP*

P2,0);

IRPP2=max(rIRPP2*P2*IRP,0);

P3RP=max(mP3*P3,0);P3H3=gP3H3*P3*H3;RPP3=max(gRPP3*RP*P3,0);

IRPP3=max(rIRPP3*P3*IRP,0);

if IRP<=0

 IRPP2=0;

IRPP3=0;

elseif IRP-IRPP2-IRPP3-max(IRP*mIRPRP,0)+RPIRP<0

if P2~=0

 IRPP2=rIRPP2*(IRP-max(IRP*mIRPRP,0)+RPIRP)/(rIRPP2+rIRPP3);

end

 if P3~=0

 IRPP3=rIRPP3*(IRP-max(IRP*mIRPRP,0)+RPIRP)/(rIRPP2+rIRPP3);

end 

end 

if P2+IRPP2+RPP2-P2RP-P2H2-P2H3-P2H1<belownoreproduction

if P2+IRPP2+RPP2-P2RP<belownoreproduction

P2RP=P2+IRPP2+RPP2;

P2H2=0;P2H3=0;P2H1=0;

else

 totP2demand=P2H2+P2H3+P2H1;

P2avail=P2+IRPP2+RPP2-P2RP;

P2H2=P2H2*P2avail/totP2demand;

P2H3=P2H3*P2avail/totP2demand;

P2H1=P2avail-(P2H2+P2H3);

end 

end 

% 

% P3 

% 

if P3+IRPP3+RPP3-P3RP-P3H3<belownoreproduction

if P3+IRPP3+RPP3-P3RP<belownoreproduction

P3RP=P3+IRPP3+RPP3;

P3H3=0;

else

 totP3demand=P3H3;

P3avail=P3+IRPP3+RPP3-P3RP;

P3H3=P3H3*P3avail/totP3demand;

end 

end 

% 

% H1 

% 

H1RP=max(mH1*H1,0);

H1HH=H1HHdemand*numHH;

if H1+P1H1+P2H1-H1RP-H1C1-H1HH<0

 if H1+P1H1+P2H1-H1RP<0

 H1RP=H1+P1H1+P2H1;

H1C1=0;H1HH=0;

else

 totH1demand=H1C1+H1HH;

H1avail=H1+P1H1+P2H1-H1RP; 

80
 H1C1=H1avail*H1C1/totH1demand;

H1HH=H1avail-H1C1;

end 

else

 if H1massdeficit<0

 H1HH=H1HH+min(H1+P1H1+P2H1-H1RP-H1C1-H1HH,-H1massdeficit);

%only what you need to make up deficit

end 

end 

% 

% H2 

% 

H2C1=gH2C1*C1*H2;H2C2=gH2C2*H2*C2;H2RP=max(mH2*H2,0);

if H2+P1H2+P2H2-H2RP-H2C1-H2C2<belownoreproduction

if H2+P1H2+P2H2-H2RP<belownoreproduction

H2RP=H2+P1H2+P2H2;

H2C1=0;H2C2=0;

else

 totH2demand=H2C1+H2C2;

H2avail=H2+P1H2+P2H2-H2RP;

H2C1=H2C1*H2avail/totH2demand;

H2C2=H2avail-H2C1;

end 

end 

% 

% H3 

% 

H3RP=max(mH3*H3,0);H3C2=gH3C2*H3*C2;

if H3+P2H3+P3H3-H3RP-H3C2<belownoreproduction

if H3+P2H3+P3H3-H3RP<belownoreproduction

H3RP=H3+P2H3+P3H3;

H3C2=0;

else

 totH3demand=H3C2;

H3avail=H3+P2H3+P3H3-H3RP;

H3C2=H3C2*H3avail/totH3demand;

end 

end 

% 

% C1 

% 

C1RP=max(mC1*C1,0);

if C1+H1C1+H2C1-C1RP<belownoreproduction

C1RP=C1+H1C1+H2C1;

end 

% 

% C2 

% 

C2RP=max(mC2*C2,0);

if C2+H2C2+H3C2-C2RP<belownoreproduction

C2RP=C2+H2C2+H3C2;

end 

81
% 

% HH 

% 

HHRP=ceil(mHH*numHH)*percapmass; 

% 

% RP 

% 

IRPRP=max(IRP*mIRPRP,0);

RPIS=min(lambda*P1IS/theta,RPISdemand);

stockRP=RP+P1RP+P2RP+P3RP+H1RP+H2RP+H3RP+C1RP+C2RP+HHRP+IRPRP;

if stockRP<0

 stockRP=0;

end 

if stockRP-(RPP1+RPP2+RPP3)-RPIRP-RPIS<=0&RPIRP==0 %this line changed

% 26 Aug 2004

RPdemand=RPP1+RPP2+RPP3+RPISdemand;

RPP1=RPP1*stockRP/RPdemand;

RPP2=RPP2*stockRP/RPdemand;

RPP3=RPP3*stockRP/RPdemand;

if RPIS~=0

 RPIS=stockRP-(RPP1+RPP2+RPP3);

else

 RPIS=0;

end 

end 

P1IS=min(theta*RPIS/lambda,P1IS); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% make checks again, to balance flows

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% P1 

% 

if P1+RPP1-P1RP-P1H2-P1H1-P1HH-P1IS<0 %if statement to deal with going

% negative

if P1+RPP1-P1RP<0

 P1RP=P1+RPP1;

P1H2=0;P1H1=0;P1HH=0;P1IS=0;

else

 totP1demand=P1H2+P1H1+P1HH+P1IS;

P1avail=P1+RPP1-P1RP;

P1H2=P1avail*P1H2/totP1demand;

P1H1=P1avail*P1H1/totP1demand;

P1HH=P1avail*P1HH/totP1demand;

P1IS=P1avail-(P1H2+P1H1+P1HH);

end 

else

 if P1massdeficit<0 %if there is an accumulated deficit between 

% demand for P1 try to

%make this up if there is extra stock

P1surplus=min(P1+RPP1-P1RP-P1H2-P1H1-P1HH-P1IS,-P1massdeficit);

%only what you need to make up deficit 

82
 P1H1=P1H1+P1surplus*P1H1massdeficit/P1massdeficit;

P1IS=P1IS+P1surplus*P1ISmassdeficit/P1massdeficit;

P1HH=P1HH+P1surplus*P1HHmassdeficit/P1massdeficit;

end

end 

% 

% P2 

% 

if IRP<=0

 IRPP2=0;

IRPP3=0;

elseif IRP-IRPP2-IRPP3-max(IRP*mIRPRP,0)+RPIRP<0

if P2~=0

 IRPP2=rIRPP2*(IRP-max(IRP*mIRPRP,0)+RPIRP)/(rIRPP2+rIRPP3);

end

 if P3~=0

 IRPP3=rIRPP3*(IRP-max(IRP*mIRPRP,0)+RPIRP)/(rIRPP2+rIRPP3);

end 

end 

if P2+IRPP2+RPP2-P2RP-P2H2-P2H3-P2H1<belownoreproduction

if P2+IRPP2+RPP2-P2RP<belownoreproduction

P2RP=P2+IRPP2+RPP2;

P2H2=0;P2H3=0;P2H1=0;

else

 totP2demand=P2H2+P2H3+P2H1;

P2avail=P2+IRPP2+RPP2-P2RP;

P2H2=P2H2*P2avail/totP2demand;

P2H3=P2H3*P2avail/totP2demand;

P2H1=P2avail-(P2H2+P2H3);

end 

end 

% 

% P3 

% 

if P3+IRPP3+RPP3-P3RP-P3H3<belownoreproduction

if P3+IRPP3+RPP3-P3RP<belownoreproduction

P3RP=P3+IRPP3+RPP3;

P3H3=0;

else

 totP3demand=P3H3;

P3avail=P3+IRPP3+RPP3-P3RP;

P3H3=P3H3*P3avail/totP3demand;

end 

end 

% 

% H1 

% 

if H1+P1H1+P2H1-H1RP-H1C1-H1HH<0

 if H1+P1H1+P2H1-H1RP<0

 H1RP=H1+P1H1+P2H1;

H1C1=0;H1HH=0;

else

 totH1demand=H1C1+H1HH; 

83
 H1avail=H1+P1H1+P2H1-H1RP;

H1C1=H1avail*H1C1/totH1demand;

H1HH=H1avail-H1C1;

end 

else

 if H1massdeficit<0

 H1HH=H1HH+min(H1+P1H1+P2H1-H1RP-H1C1-H1HH,-H1massdeficit);

%only what you need to make up deficit

end 

end 

% 

% H2 

% 

if H2+P1H2+P2H2-H2RP-H2C1-H2C2<belownoreproduction

if H2+P1H2+P2H2-H2RP<belownoreproduction

H2RP=H2+P1H2+P2H2;

H2C1=0;H2C2=0;

else

 totH2demand=H2C1+H2C2;

H2avail=H2+P1H2+P2H2-H2RP;

H2C1=H2C1*H2avail/totH2demand;

H2C2=H2avail-H2C1;

end 

end 

% 

% H3 

% 

if H3+P2H3+P3H3-H3RP-H3C2<belownoreproduction

if H3+P2H3+P3H3-H3RP<belownoreproduction

H3RP=H3+P2H3+P3H3;

H3C2=0;

else

 totH3demand=H3C2;

H3avail=H3+P2H3+P3H3-H3RP;

H3C2=H3C2*H3avail/totH3demand;

end 

end 

% 

% C1 

% 

if C1+H1C1+H2C1-C1RP<belownoreproduction

C1RP=C1+H1C1+H2C1;

end 

% 

% C2 

% 

if C2+H2C2+H3C2-C2RP<belownoreproduction

C2RP=C2+H2C2+H3C2;

end 

84
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

ISHHflow=max(0,(theta+lambda)*ISHHdemand*numHH);

ISIRP=ISHHflow;

if ISmass+P1IS+RPIS-ISIRP<=0

 ISIRP=ISmass+P1IS+RPIS;

else

 if ISmassdeficit<0&numHH>=2 %if there is an accumulated deficit 

% between demand for IS by the HH and IS

%supplied, try to make this up if there is extra stock

ISIRP=ISIRP+min(ISmass+P1IS+RPIS-ISIRP,-ISmassdeficit);%only

% what you need to make up deficit

end 

end 

if (P1HH+H1HH+ISIRP)==0

weightedprice=0;

percapbirths=0;

elseif (pP1*P1HH+pH1*H1HH+pIS*ISIRP)==0

weightedprice=0;

percapbirths=0;

else 

%weightedprice=(pP1*P1HHdemand+pH1*H1HHdemand+pIS*ISHHflow)/(P1HHdemand

% +H1HHdemand+ISHHflow);

weightedprice=(pP1*P1HH+pH1*H1HH+pIS*ISIRP)/(P1HH+H1HH+ISIRP);

%percapbirths=max(etaa-

% etab*W/weightedprice+etac*(W/weightedprice)^2,0);

percapbirths=max(etaa-etab*sqrt(W/weightedprice),0);

% if etab<=2*etac*(W/weightedprice)|etaa<=etac*(W/weightedprice)^2

% error('pricesettingmodel:growtherror','Percapbirths function

% parameters give bad function form')

% end 

end 

% 

%%%%% changed 26 AUgust 2004

% 

nextP1 = P1+RPP1-P1RP-P1H2-P1H1-P1HH-P1IS; 

if P1==0

 P1H1demand=0;P1ISdemand=0;P1HHdemand=0;

end 

nextP1H1massdeficit=P1H1massdeficit+P1H1-P1H1demand; 

nextP1ISmassdeficit=P1ISmassdeficit+P1IS-P1ISdemand; 

nextP1HHmassdeficit=P1HHmassdeficit+P1HH-P1HHdemand*numHH; 

85
nextP2 = P2+IRPP2+RPP2-P2RP-P2H2-P2H3-P2H1;

nextP3 = P3+IRPP3+RPP3-P3RP-P3H3;

nextH1 = H1+P1H1+P2H1-H1RP-H1C1-H1HH;

if H1==0

 H1HHdemand=0;

end 

nextH1massdeficit=H1massdeficit+H1HH-H1HHdemand*numHH; 

nextH2 = H2+P1H2+P2H2-H2RP-H2C1-H2C2;

nextH3 = H3+P2H3+P3H3-H3RP-H3C2;

nextC1=C1+H1C1+H2C1-C1RP;

nextC2=C2+H2C2+H3C2-C2RP;

nextHH = HH+P1HH+H1HH-HHRP;

nextISmass=ISmass+P1IS+RPIS-ISIRP; %keep track of actual mass in IS

nextISmassdeficit = ISmassdeficit+ISIRP-ISHHflow;

%keep track of deficit in IS, what is supplied minus the demand

nextIRP = IRP-IRPP2-IRPP3+RPIRP+ISIRP-IRPRP; %this line changed 2 Aug

% 2004

nextRP = stockRP-(RPP1+RPP2+RPP3)-RPIRP-RPIS;

nextnumHH=max(numHH+ceil(percapbirths*numHH)-ceil(mHH*numHH)-

ceil(numHH*phi*(percapmass-idealpercapmass)^2),1);

nextpercapmass=nextHH/nextnumHH;

sys =

[nextP1;nextP2;nextP3;nextH1;nextH2;nextH3;nextC1;nextC2;nextHH;nextISm

ass;nextRP;nextIRP;... 

nextP1H1massdeficit;nextP1ISmassdeficit;nextP1HHmassdeficit;nextH1massd

eficit;nextISmassdeficit;...

nextnumHH;nextpercapmass]; 

% end mdlUpdate 

%

%======================================================================

% mdlOutputs

% Return the block outputs.

%======================================================================

%

86
function sys=mdlOutputs(t,y,u)

global RPP1 P1H2 P1IS P1H1 P1HH RPIS ISIRP P2H1 H1C1 H1HH pP1 pH1 pIS W

ISHHflow

sys =

[y;RPP1;P1H2;P1IS;P1H1;P1HH;RPIS;ISIRP;P2H1;H1C1;H1HH;pP1;pH1;pIS;W];

% end mdlOutputs

%

%======================================================================

% mdlGetTimeOfNextVarHit

% Return the time of the next hit for this block. Note that the result 

% is absolute time. Note that this function is only used when you 

% specify a variable discrete-time sample time [-2 0] in the sample 

% time array in

% mdlInitializeSizes.

%======================================================================

%

function sys=mdlGetTimeOfNextVarHit(t,y,u)

sampleTime = 1; % Example, set the next hit to be one second later.

sys = t + sampleTime;

% end mdlGetTimeOfNextVarHit

%

%======================================================================

% mdlTerminate

% Perform any end of simulation tasks.

%======================================================================

%

function sys=mdlTerminate(t,y,u)

sys = [];

% end mdlTerminate

87
References 

Barbier, E. B., and A. Markandya, “The conditions for achieving environmentally 

sustainable development.” In The Economics of Environment and Development (E. B. 

Barbier, ed), Cheltenham, UK: Edward Elgar Publishing, 1998, 22-42. 

Boumans, R., Costanza, R., Farley, J., Villa, F. and Wilson, M. "Modeling the dynamics 

of the integrated earth system and the value of global ecosystem services using the 

GUMBO model," Special Issue: The Dynamics and Value of Ecosystem Services: 

Integrating Economic and Ecological Perspectives, Ecological Economics, 41, (2002) 

529-560.. 

Bird, P. J. W. N. The transferability and depletability of externalities. Journal of 

Environmental Economics and Management 14:54-57 (1987). 

Bond, E. M. and J. M. Chase. Biodiversity and ecosystem functioning at local and 

regional spatial scales. Ecology Letters 5:467-470 (2002). 

Brock, W., and A. Xepapadeas. Optimal ecosystem management when species compete 

for limiting resources. Journal of Environmental Economics and Management 44(2):189­

220 (2002). 

Brock, W.A., and D. Starrett. Managing systems with non-convex positive feedback. 

Environmental and Resource Economics 26:575-602 (2003). 

Cabezas, H., C. W. Pawlowski, A. L. Mayer, and N. T. Hoagland. Sustainability: 

ecological, social, economic, technological, and systems perspectives. Clean 

Technologies and Environmental Policy 5:167-180 (2003). 

Cabezas, H., C. W. Pawlowski, A. L. Mayer, and N. T. Hoagland. Sustainable systems 

theory: ecological and other aspects. Journal of Cleaner Production 13(5):455-467 

(2005). 

Cabezas, H., C.W. Pawlowski, H. Wm. Whitmore, and A.L. Mayer. On the sustainability 

of integrated model systems with industrial, ecological, and macroeconomic components. 

Resources, Conservation, and Recycling, (2006), in press. 

Carpenter, S. R., and W. A. Brock. Spatial complexity, resilience, and policy diversity: 

fishing on lake-rich landscapes. Ecology and Society 9(1):8 (2004). 

Carpenter, S., W. Brock, and P. Hanson. Ecological and social dynamics in simple 

models of ecosystem management. Conservation Ecology 3(2) (1999). 

Carpenter, S., B. Walker, J. M. Anderies, and N. Abel. From metaphor to measurement: 

resilience of what to what? Ecosystems 4:765-781 (2001). 

88
Conrad, Klaus. Computable general equilibrium models in environmental and resource 

economics. In The International Yearbook of Environmental and Resource Economics 

2002/2003 (T. Tietenberg and H. Folmer, eds.), Edward Elgar, Northampton, MA. 2002, 

Costanza, R., J. Cumberland, H. Daly, R. Goodland, and R. Norgaard. An Introduction to 

Ecological Economics. Boca Raton, FL: St. Lucie Press, 1997. 

Costanza, R., R. D'Arge, R. de Groot, S. Farber, M. Grasso, B. Hannon, K. Limburg, S. 

Naeem, R. V. O'Neill, J. Paruelo, R. G. Raskin, P. Sutton, and M. van den Belt. The 

value of the world’s ecosystem services and natural capital. Nature 387:253-260 (1997). 

Costanza, R., A. Voinov, R. Boumans, T. Maxwell, F. Villa, L. Wainger, and H. Voinov. 

Integrated ecological economic modeling of the Patuxent River watershed, Maryland. 

Ecological Monographs 72(2):203-231 (2002). 

Daily, G. C. Nature’s Services: Societal Dependence on Natural Ecosystems. Washington 

D.C.: Island Press, 1997. 

Dasgupta, P., S. Levin, and J. Lubchenco. Economic pathways to ecological 

sustainability. BioScience 50(4):339-345 (2000). 

Diamond, J. Guns, Germs, and Steel: The Fates of Human Societies. New York/London: 

W.W. Norton and Company, 1999. 

Doak, D. F., and L. S. Mills. A useful role for theory in conservation. Ecology 75:615­

626 (1994). 

Ehrlich, P. R. The scale of the human enterprise and biodiversity loss. In Extinction 

Rates (J. H. Lawton, and R. M. May, eds.). Oxford, UK: Oxford University Press, 1995, 

215-226. 

Farber, S. C., R. Costanza, and M. A. Wilson. Economic and ecological concepts for 

valuing ecosystem services. Ecological Economics 41:375-392 (2002). 

Fath, B. D., H. Cabezas, and C. W. Pawlowski. Regime changes in ecological systems: 

an Information Theory approach. Journal of Theoretical Biology 222(4):517-530 (2003). 

Ferenc, J., S. Istvan, V. Gabor. Species positions and extinction dynamics in simple food 

webs. Journal of Theoretical Biology 215:441-448 (2002). 

Freeman, A. M. Depletable externalities and Pigouvian taxation. Journal of 

Environmental Economics and Management 11:173-179 (1984). 

Gatto, M. Sustainability: is it a well defined concept? Ecological Applications 5(4):1181­

1183 (1995). 

89
Goodland, R., and H. Daly. Environmental sustainability: universal and non-negotiable. 

Ecological Applications 6(4):1002-1017 (1996). 

Gunderson, L. H. Ecological resilience - in theory and application. Annual Review of 

Ecology and Systematics 31:425-439 (2000). 

Haberl, H., M. Wackernagel, F. Krausmann, K. H. Erb, and C. Monfreda. Ecological 

footprints and the human appropriation of net primary production: a comparison. Land 

Use Policy 21:279-288 (2004). 

Hardin, G. The tragedy of the commons. Science 162(1968):1243-1248 (1968). 

Heal, G. Nature and the Marketplace: Capturing the Value of Ecosystem Services. 

Washington, D.C.: Island Press, 2000. 

Heller, M. A. The three faces of private property, 79 Oregon Law Review 417 (2000). 

Heller, M. A. The tragedy of the anticommons: property in the transition from Marx to 

markets, 111 Harvard Law Review 621 (1998). 

Holling, C. S. Resilience and stability of ecological systems. Annual Review of Ecology 

and Systematics 4:1-23 (1973). 

Holling, C. S. Surprise for science, resilience for ecosystems, and incentives for people. 

Ecological Applications 6(3):733-735 (1996). 

Hooper, D. U., F. S. Chapin, J. J. Ewel, A. Hector, P. Inchausti, S. Lavorel, H. H. 

Lawton, D. M. Lodge, M. Loreau, S. Naeem, B. Schmid, H. Setala, A. J. Symstad, J. 

Vandermeer, and D. A. Wardle. Effects of biodiversity on ecosystem functioning: A 

consensus of current knowledge. Ecological Monographs 75(1):3-35 (2005). 

Imhoff, M. L., L. Bounoua, T. Ricketts, C. Loucks, R. Harriss, and W. T. Lawrence. 

Global patterns in human consumption of net primary productivity. Nature 429:870-873 

(2004). 

Kearns, C. A. Pollinators, flowering plants, and conservation biology. BioScience 47:297­

307 (1997). 

Krutilla, J. Conservation reconsidered. American Economic Review 57:777-786 (1967). 

Levin, S. A., J. Dushoff, and J. E. Keymer. Community assembly and the emergence of 

ecosystem pattern. Scientia Marina 65 (Supp. 2):171-179 (2001). 

Ludwig, D., S. Carpenter, and W. Brock. Optimal phosphorus loading for a potentially 

eutrophic lake. Ecological Applications 13(4):1135-1152 (2003). 

90
MacArthur, R. H., and E. O. Wilson. The Theory of Island Biogeography. Monographs 

in Population Biology. Princeton, NJ: Princeton University Press, 1967. 

Mayer, A. L., and M. Rietkerk. The dynamic regime concept for ecosystem management 

and restoration. BioScience 54(11):1013-1020 (2004). 

McMichael, A. J., C. D. Butler, and C. Folke. New visions for addressing sustainability. 

Science 302:1919-1920 (2003). 

Millennium Ecosystem Assessment Synthesis Report. Washington, D.C. World 

Resources Institute, 2005. 

O’Neill, R. V., and J. R. Kahn. Homo economus as a keystone species. BioScience 

50(4):333-337 (2000). 

Pauly, D., V. Christensen, S. Guénette, T. J. Pitcher, U. R. Sumaila, C. J. Walters, R. 

Watson, and D. Zeller. Towards sustainability in world fisheries. Nature 418:689-695 

(2002). 

Pimm, S. L., and R. A. Askins. “Forest losses predict bird extinctions in eastern North 

America.” In Proceedings of the National Academy of Sciences, USA 92:9343-9347, 

1995. 

Portela, R., and I. Rademacher. A dynamic model of patterns of deforestation and their 

effect on the ability of the Brazilian Amazonia to provide ecosystem services. Ecological 

Modelling 143:115-146 (2001). 

Postel, S., and S. Carpenter. “Freshwater ecosystem services.” In Nature’s Services: 

Societal Dependence on Natural Ecosystems. Washington D.C.: Island Press, 1997, 195­

214. 

Reed, W. J. Optimal harvesting models in forest management – a survey. Natural 

Resource Modeling 1:55-80 (1986). 

Rees, W. E. An ecological economics perspective on sustainability and prospects for 

ending poverty. Population and Environment 24(1):15-45 (2002). 

Sakai, A. K., F. W. Allendorf, and J. S. Holt. The population biology of invasive species. 

Annual Review of Ecology and Systematics 32:305-332 (2001). 

Samuelson, P. The pure theory of public expenditure. Review of Economics and Statistics 

36:387-389 (1954). 

Scheffer, M., S. Carpenter, J. A. Foley, C. Folke, and B. Walker. Catastrophic shifts in 

ecosystems. Nature 413(6856):591-596 (2001). 

91
Settle, C., T. D. Crocker, and J. F. Shogren. On the joint determination of biological and 

ecological systems. Ecological Economics 42:301-311 (2002). 

Slade, M. E. Trends in natural resource commodity prices: An analysis of the time 

domain. Journal of Environmental Economics and Management 9:122-137 (1982). 

Tilman, D. The ecological consequences of changes in biodiversity: A search for general 

principles. Ecology 80:1455-1474 (1999). 

United Nations Development Programme. URL: 

www.unfpa.org/swp/2001/english/ch03.html (1998). 

United States Census Bureau. URL: www.census.gov/ipc/www/worldpop.html 

(2005). 

van den Bergh, J. J. M. ”A multisectoral growth model with materialsflows and 

economic-ecological interactions.” In Ecological Economics and Sustainable 

Development, Northampton, MA: Edward Elgar Publishing, 1996. 

Voinov, A., R. Costanza, L. Wainger, R. Boumans, F. Villa, T. Maxell, and H. Voinov. 

Patuxent landscape model: integrated ecological economic modeling of a watershed. 

Environmental Modelling and Software 14:473-491 (2002). 

Wardle, D. A., K. I. Bonner, and G. M. Barker. Stability of ecosystem properties in 

response to above-ground functional group richness and composition. Oikos 89:11-23 

(2000). 

World Commission on Environment and Development. 1987. Our Common Future. 

Oxford: Oxford University Press, 1987. 

Yandle, B. Grasping for the heavens: 3-D property rights and the global commons. 10 

Duke Envirmental Law and Policy Forum 13 (1999). 

Yandle, B., and A. P. Morriss. The technologies of property rights: choice among 

alternative solutions to tragedies of the commons. 28 Ecology Law Quarterly 123 (2001). 

92