Continuous Dynamic System Simulation (STB-8: ssi3) ------------------------------------ The syntax of the ^simula^ command is: . ^simula^, [^dt(^#^)^] [^tstart(^#^)^] [^tspan(^#^)^] [^exovar(^filename^)^] [^parv(^pname1=# pname2=# ....^)^] [^ival(^svar1=# svar2=# ....^)^] [graph_options] The options are: ^dt(^#^)^ specifies the integration time step (time units). The default is 1. ^tstart(^#^)^ specifies the time at which the simulation start. The default value is 0. ^tspan(^#^)^ specifies the simulation length. The default value is 30. ^exovar^ requests a file that will be merged to the model file in memory. This file contains the exogenous variables used in the model. This option is needed when an exogenous variable statement is declared in _model. ^parv^ allows the modification of parameter values, overriding the values specified in the _model statements. ^pname^ is a parameter name. ^ival^ allows the modification of the initial values of the state variables, overriding the values specified in the _model statements. ^sname^ is a name of a state variable. graph_options all the graphic option allowed for the two-ways scatterplot. For example: ^ylog^ could be useful when the expected increase is explosive, ^saving^ can be useful for saving the resulted graphic, and so on. filename name of a *.dta file. pname# name of parameter (Stata variable). svar# name of a state variable in the model (Stata variable) Notes ------- ^simula^ requires a model file already in memory; on request, it loads and merges and exogenous variable. The model can also be inserted directly with the command: ^. input str## _model^. It performs the numerical integration of the model equations by the Euler (rectangular) method. The command gives the change in time of the system conditions (state variables) starting from the initial values and with the provided parameter values. At the end of simulation it draws a graph of the state variables versus time and, for all the time step required, the state, rate and auxiliary variable values, the exogenous variables and the _time variable will be available in memory. Wishing to modify the relations or parameters from the Stata command level, is possible to use the ^replace _model=exp in #^ command to modify the model list. - Time solution step (dt) To ensure a stable solution, a good rule is to keep dt between 1/2 and 1/10 of the smallest time coefficient in the model (Richardson and Pugh, 1981). Frequently, it is suggests a ^dt^ value of 1/4-1/5 of the smallest time coefficient (Ferrari, 1978). Values higher than 1/2 determine computational inaccuracies while values lower than 1/10 require much of computer time. Moreover, a too short dt value can also increase the round-off error because the rate values added to the state variables becomes very small. Anyhow, is useful, when performing a simulation, to try with a smaller dt value and to check the amount of error. The solution is correct if a small change in the time step does not affect the simulation results. The choices for the time step and the starting time must allow the merging of the exogenous variables; the _time variable created by ^simula^ must have the values present in the _time variable of the exovar file (not only those, however). If the two _time variable values do not match, the exogenous variable's values are not merged. For example: let the _time values in the exovar file be 0, 1, 2, 3 and 4 and let choice dt=0.4, tspan=8 and tstart=1. The values for the _time variable created by ^simula^ will be 1, 1.4, 1.8, 2.2, 2.6, 3.0, 3.4 and 3.8. As effect, the program merges only the exovar values corresponding to time 1 and 3. - Exogenous variables When the exogenous variable's data are available with a lower frequency than the specifies time solution step, the command generates the missing data. The initial and final values are generated by repetition of the nearest value, the intermediate ones by linear interpolation. - Parameter and initial values It is possible to override the parameter and initial values coded in the model file by specifying them in the command options (parv and ival). Indeed, while is a good practice to indicate such values in the model file, using the parv and ival options provides a faster method to analyze the model sensitivity and the effects of changing the initial conditions. - Delays The delay functions apply to all the state variables of the equation. If in a rate equation there are state variables with different delays time, it is better to decompose the rate equation by using auxiliary equations. Note that the effect of IDL and MDL is the same until the delay time (dly) is a constant. If dly is changing, the retarded information is lost while the retarded material is conserved. Only the rate equations may contain the delay functions. - Auxiliary equations If an auxiliary variable (ex: avar1) depends on another auxiliary equation (ex: avar2) the equation for avar2 must precede the avar1 equation in the model listing. Example ------- For using ^simula^ is necessary to have a model coded in the string variable _model, present in memory. With the example on prey-predator interactions in a ecosystem, we can load the file ^volmod1.dct^ and use ^simula.ado^: ^. infile using volmod1 dictionary { * file "volmod1.dct" * example model for "simula" STATA command (4/11/91) * This is a kind of a VOLTERRA model - without delays * Ref. J. Maynard Smith, 1974. Models in ecology. Cambridge University Press. str70 _model "Model statements" } (25 observations read) If, after a reexamines of experimental results (or if we simply want to see the effects of some parameter modifications), we decide to modify the initial or parameter values, it is not necessary to modify _model but we can change them at the command level with the options ^ival^ and ^parv^. Suppose that the initial values for X and Y are 800 and 400 respectively and that Kbx=0.07 and Kby=0.003. The command will be: ^. simula,dt(0.25) tspan(10) tstart(1991) exo(climate) par(Kbx=0.07 Kby=0.003) ^> ival(X=800 Y=400) bord c(ll) Parsing: Line: 1 Line: 2 Line: 3 . . . Line: 23 Line: 24 Line: 25 Simulation: Time=1991 Time=1991.25 Time=1991.5 . . . Time=2000.25 Time=2000.5 Time=2000.75 A graph of the evolution in time of the state variables is shown. Model development ----------------- The modeling process consists in the following steps: 1) identification of the relevant system variables for the representation of the system conditions (state variables); 2) identification of the fluxes (rate variables) which determine the ^speed^ of the state variable variations. ^State equations^ link the rate variables to the state variables; 3) development of the relations for the rates. Rates are computed as function of state and exogenous variables (not other rates). For a clearer analysis of the system, the rate equations can be also function of auxiliary variables. The third step is the most important for the modeling process. In fact, the evolution of the state variables derives from the integration of the rate equations, starting from certain ^initial conditions^. Example: Let we consider a classical Volterra model on prey-predator interactions in an ecological system (see Maynard, 1974 for details) and assume that the prey is a herbivorous. Two state variables describe the system conditions: the prey density in the ecosystem (X, as number of individuals per unit area) and its predator density (Y, as number of individuals per unit area). For the model development we made the following assumptions: 1) The birth rate of the prey, without predator and without limiting factors, is a constant fraction of the prey density. 2) The environmental limits determine a maximum value for the prey density that is not possible to exceed. This limit is generally indicated as ^carrying capacity^ (Cc). 3) The predation rate is a fixed fraction of all the possible meets between prey and predator (represented by the product X*Y). 4) The death rate of predator is a constant fraction of the predator density. 5) The birth rate of predator is a fraction of Y density but the predator birth rate tends to zero when the X density approach to the Y density. 6) In addition, to show the possibility of use the exogenous variables, we make a further assumption: if the prey is a herbivorous, the carrying capacity of the environment (Cc) depends on the annual grass production of the range. Cc, in turn, could be dependent on the annual ratio (total rainfalls/average temperature) by the equation: Cc=500+K1*rain/temp. ^rain^ and ^temp^ are exogenous variables. The differential equations describing the change in time of X and Y are: 1) dX/dt = Kbx*X * (Cc-X)/Cc - Kmx*(X*Y) variation X effect of death rate rate of birth carrying for X rate capacity predation 2) dY/dt = -Kmy*Y + Kby*Y*(X-Y) variation Y effect of predation rate of death on predator birth Y rate where 1/Kbx, 1/Kmx, 1/Kmy and 1/Kby are the time coefficients of the rates??. Let the initial values for prey and predator density in the area be 1000 and 100 individuals, respectively. Kbx is the birth coefficient of the prey population, expressed as fraction of X birth per year, and fixed to 0.05. The birth coefficient for predator is Kby=0.006 and represents the efficiency of the conversion prey-predator. The mortality coefficients for prey and predator are, respectively, Kmx=0.001 and Kmy=0.05. The coefficient for grass production (K1) equals to 80. Models coding ------------- The model equations have to be written in a form understandable by the ^simula^ command. ^simula^ uses two kinds of Stata file: the ^model file^ and the ^exogenous variable file^. The models are to be saved in the ^_model^ variable of the model file; _model is a string variable and can have up to 80 characters. There are six kinds of statements recognized by ^simula^: state equations (S), auxiliary equations (A), rate equations (R), exogenous variable's declaration (E), initial value's assignment to the state variables (I) and parameter value's assignment (P). The above symbols (S, A, R, E, I and P), placed as first character of each line identifies the statement type. ^simula^ parsed the statements in a way requiring the following rules: 1) Line structure Each statement is composed by three parts: a) a statement identifier (the capital letters S, A, R, E, I or P); b) the statement itself (equation, declaration, assignment); blanks inside are not allowed; c) a statement label (is optional and must be separated by a blank, at least, from the statement). 2) Comment lines Any group of characters different from S, A, R, E, I and P denotes a "comment line". 3) Lines order The model list can contain the statements in any order The only exception is when an auxiliary variable (aux1) depends on another auxiliary variable (aux2): in this case the equation for aux2 must precede the aux1 equation. 4) All the statement kinds require only one expression/declaration/assigment per line. For example, it's not possible to write the lines: E rain temp, I X=100 Y=50, P K1=5 K2=7; they should be rewritten as E rain, E temp, P K1=5, P K2=7. 5) Constant values The S, A, R, I and P statements allow constant values. Arguments for the special functions are to be parameters or auxiliary variable (ex: the delay time) and not constant. 6) Special functions Right now ^simula^ has three special functions: IDL (information delay), MDL (material delay) and SWTC (switch a rate on/off). Functions can be only specified in the rate equations; arguments must be parameters or auxiliary variables. 7) Blank lines Between two statements is possible to insert only one blank line. This insertion could be used to make a clearer reading of the model. 8) Parameters and initial values It is not necessary to declare in _model all the parameter and initial values, but this is a good practice. It is possible to assign those values as option in the ^simula^ command. Of course, all parameter and initial values require its declaring; this can be done or in _model or as option. 9) Statement syntax The syntax, for each kind of statement, is the following: S svar1=svar1+dt*(rvar1+rvar2+rvar3+...) label A avar=exp label R rvar=exp label E evar label I svar=# label P pname=# label where: svar state variable name avar auxiliary variable name rvar rate variable name evar exogenous variable name declared pname parameter name # numeric constant dt time step for the equation integration. It is inserted only for identifying a state equation but, really, any character set differing by ^+^ or ^*^ is allowed. exp any algebraic expression containing #, pname, evar, avar and svar. In the expressions, all the Stata functions can also be used. For example is possible to use the random generator uniform(). label a descriptor of the statement (optional). 10) Special functions ^simula^ recognizes some special functions that, indicated in the rate equations, perform some specific tasks. These can to be indicated only in the rate equations: moreover, each rate equation allows only one function. The syntax of the functions is: R rvar=fnc[exp,arg,arg,...] where ^fnc^ is a generic function, ^exp^ an algebraic expression and ^arg^ are the arguments needed by the function itself. The available functions are, at present: a) IDL information delay function. The syntax is: ^R rvar=IDL[exp,dly]^ where ^dly^ is the delay time. The effect on rvar is the follow: rvar(t)=exp(t-dly) if t>=dly rvar(t)=0 if t>dly t is the time index. b) MDL material delay function; this function has the same effect of a further level inserted in the flux, but not available to the user. The syntax is: ^R rvar=MDL[exp,dly]^,where ^dly^ is the delay time. The effect on rvar is the follow: rvar(t)=rvar(t)+exp(t-dly) if t>=dly rvar(t)=0 if t>dly ^dly^ is expressed in time units. It must be a parameter or an auxiliary variable name (not a constant or expression). c) SWTC switch to rate on/off. Equals the rate to zero if a ^control^ ^variable^ (cvar) is below a low threshold value (von) or reach an upper threshold value (voff). The syntax is: ^R rvar=SWTC[exp,cvar,von,voff]^ and the results are: rvar=0 if cvar<=von rvar=exp if cvar>von and cvar=voff This function performs different tasks: 1) If a physical constraint limits a state variable and the inflow rate is independent from the level: von is set to 0 and voff to the level maximum capacity; cvar is the level itself. 2) When a process start and/or stops at some defined times: the control variable is _time, ^von^ and ^voff^ are the starting and ending times. 3) Gives a pulse of material at a specific time only. 4) Other processes conditionally associated to other variables. Example: Under, the *.DCT file containing the Volterra example model for the ^simula^ command, is reported. dictionary { * file "volmod1.dct" * example model for "simula" STATA command (4/11/91) * This is a kind of a VOLTERRA model - without delay * Ref. J. Maynard Smith, 1974. Models in ecology. Cambridge University Press. str70 _model "Model statements" } "* Prey Predator Volterra model - without delay " "* STATE equations ----------- " "S X=X+dt*(Bx-Dx) Prey population (n) " "S Y=Y+dt*(-Dy+By) Predator population (n) " "* RATE equations ------------ " "R Bx=Kbx*(Cc-X) Birth rate of prey (n/year) " "R Dx=Kmx*X*Y Death rate of prey (n/year) " "R Dy=Kmy*Y Death rate of predator (n/year) " "R By=Kby*Y*(X-Y) Birth rate of predator (n/year) " "* AUXILIARY equations ------- " "A Cc=500+K1*rain/temp Environm. carrying capacity (n) " "* e.g. related to the ecosystem grass production " "E temp Aver. annual temperature (C) " "E rain Aver. annual rainfall (mm) " "I X=1000 Prey initial density " "I Y=100 Predator initial density " "P Kbx=0.05 Fractional birth rate of X " "P Kby=0.006 Birth coefficient for Y " "P Kmx=0.001 Predation coefficient " "P Kmy=0.05 Fractional death rate for Y " "P K1=80 Effect of the rain/temp ratio " Exogenous variable file creation -------------------------------- Exogenous variables required by the model are to be saved in an exogenous variable file, also containing a time variable named ^_time^. Before saving this file is necessary to sort it by _time because, ^simula^ merges the exogenous variable file to the model by the _time variable. Example: The exogenous variable file created contains the rainfall and temperature values for ten years. We can assume that the values came from a stochastic weather model; they are inserted in the following *.DTA file: ^. d^ Contains data from d:\climate.dta Obs: 10 (max= 14230) EXO-VARS for the Volterra model Vars: 3 (max= 99) 1. _time float %9.0g Year 2. temp float %9.0g Annual average temper. (C) 3. rain float %9.0g Total annual rainfalls (mm) Sorted by: _time ^. list^ _time temp rain 1. 1991 18 1550 2. 1992 17 1300 3. 1993 15 1800 4. 1994 21 1950 5. 1995 15 1200 6. 1996 22 1100 7. 1997 19 1500 8. 1998 20 1400 9. 1999 15 1700 10. 2000 16 1600 Author ------ Francesco Danuso, Universita degli Studi di Udine, Istituto di Produzione Vegetale, Udine, Italy. Fax, 011-39-432-558603