program define simula
version 2.1
*   CONTINUOS DYNAMIC SISTEM SIMULATION WITH STATA (v.2)
*          (Francesco Danuso, 18/2/1992)

*  (Main ref: Principles of Systems, J.W. Forrester, 1968, M.I.T. Press)
*
* sintax: . simula, [dt(#)] [tspan(#)] [tstart(#] [exovar(filename)]
*                   [parv(pname1=# pname2=# ...)] [ival(svar1=# svar2=# ...)]
*                   [graph_options]

* special function allowed: IDL, MDL, SWTC.

mac def ck "dt(real 1) TSPan(real 30) TSTart(real 0)  *"
mac def _options "EXOvar(string) PARv(string) IVal(string) %ck"

parse "%_*"

* The exogenous variable file is loaded if requested.
* NOTE:
* - the model file must be already in  memory when 'simula' is called.
* - the model string variable must be named '_model'
* - the file with exovar will be loaded if the option 'exovar' is present.
* - the file exovar must contain the exogenous variable and a time variable
*   named '_time'.

capture confirm string var _model
if _rc {
   disp in red " MODEL NOT FOUND (load or input a model)"
   exit
}


quietly {
mac def tend=%_tstart+%_tspan
mac def nobs=0.5+%_tspan/%_dt
mac def nobs=int(%nobs)
* drop if _model==""
if %nobs<_N {mac def nobs=_N}

set obs %nobs

capture drop _time
capture drop bufs
capture drop bufr
capture drop zero
gen _time=(_n-1)*%_dt+%_tstart  in 1/%nobs       /* time variable */
gen double bufs=.
gen double bufr=0
gen zero=0
lab vari _time "Time"

sort _time
}

set more 1

* ======== Parsing =====================================
disp _n(1) "Parsing: "
mac def svar ""     /*   list of state variables    */
mac def avar ""     /*   list of auxiliary variables */
mac def rvar ""     /*   list of rate variables  */
mac def evar ""     /*   list of exogenous variable  */

mac def k=1
while _model[%k]!="" | _model[%k+1]!="" {
disp in gre "Line: %k"
mac def w=_model[%k]
parse "%w", parse(" ")
mac def keq "%_1"                /* kind of equation: S/R/A/I/E/P  */
mac def eq "%_2"
* -- re-builds the statement label --
mac shift
mac shift
mac def lbl ""
while "%_1"!="" {
   mac def lbl="%lbl %_1"
   mac shift
}
parse "%w", parse(" ")

   * ------- parses STATE equations --------------------------------
   if "%keq"=="S" {
       parse "%eq",parse("=*")        /*  ex: a=a+dt*(b+c-d-f)             */
       mac def svar "%svar %_1"       /*  state variable list  */
       mac def seq%_1="%_5"         /*  %_2 =, %_5 the variation rate    */
       capture drop %_1
       qui gen double %_1=.
       label vari %_1 "%lbl"
   }

   * ------- parses AUXILIARY equations ------------------
   if "%keq"=="A" {
       parse "%eq",parse("=")            /*   ex: a=10+X/Y   */
       mac def avar "%avar %_1"
       mac def aeq%_1="%_3"
       capture drop %_1
       qui gen double %_1=.
       label vari %_1 "%lbl"
   }


   * --------- parses RATE equations ----------------------
   if "%keq"=="R" {
       parse "%eq",parse("=[,]")     /*   ex: Bx=MDL[1/10*(60-peso),dly]   */
       mac def rvar="%rvar %_1"
       mac def func "%_3"
       mac def fcod%_1=0       /* func. code 0=no, 1=info. del,2=mater. delay */
       mac def req%_1="(%_3)"  /*      3=switch on, 4=switch off.             */

       if "%func"=="IDL" {                /*  information delay required */
          mac def fcod%_1=1
          mac def req%_1="(%_5)"
          mac def idl%_1="%_7"
       }
       if "%func"=="MDL" {                /*  material delay required */
          mac def fcod%_1=2
          mac def req%_1="(%_5)"          /*  %_7 holds the delay time var. */
          mac def mdl%_1="%_7"
       }
       if "%func"=="SWTC" {             /*  switch-on/off rate */
          mac def fcod%_1=3             /*  ex R rvar=SWTC[exp,cvar,von,voff] */
          mac def req%_1="(%_5)"        /*  expression  */
          mac def cvar%_1="%_7"         /*  control variable */
          mac def von%_1="%_9"          /*  on threshold value */
          mac def voff%_1="%_11"        /*  off threshold value */
       }

       capture drop %_1
       qui gen double %_1=0
       label vari %_1 "%lbl"
   }


   * --------- parses INITIAL VALUES ------------
   if "%keq"=="I" {
       qui replace %_2 in 1    /* assigns init. values, ex: ivpeso=10 */
   }

   * ----------- parses PARAMETERS -----------------------
   if "%keq"=="P" {
       qui capture  gen  %_2      /* generates a parameter variable */
       if _rc!=0 {qui replace %_2}
   }

   * ----------- parses EXOGENOUS VARIABLES ----------------
   if "%keq"=="E" {
       mac def evar="%evar %eq"     /* list of exovar   */
   }

mac def k=%k+1            /*  index of a new model statement  */
}


* ========= Clears, merges, checks and fills the holes in EXOVARS =========

if "%evar"!="" {

    disp _n(1) "Merging file: %_exovar" _n(1)

    parse "%evar", parse(" ")
    while "%_1"!="" {
       capture drop %_1              /*  drops exovars to merge new values */
       mac shift                     /*  allowing the re-run of simulation */
    }

    if "%_exovar"!="" {
       merge _time using %_exovar     /* links exovars */
       drop _merge
    }

    parse "%evar", parse(" ")           /* checks */
    while "%_1"!="" {
         capture confirm var %_1
         if _rc {
            disp in red "The exovar '%_1' is not present."
            disp in red "Check the EXOGENOUS VAR. FILE in 'exovar' option."
            exit
         }

         * ---- fills, with linear interpolation ---
         mac def i=1
         while %i<=_N {
             if %_1[%i]==. {
                mac def t1=_time[%i-1]
                if %t1==. {mac def t1=_time[1]}
                mac def v1=%_1[%i-1]
                mac def j=%i
                while %_1[%j]==. & %j<_N {
                    mac def j=%j+1
                }
                mac def t2=_time[%j]
                mac def v2=%_1[%j]
                if %v1==. {mac def v1=%v2}
                if %v2==. {mac def v2=%v1}
                mac def Slope=(%v2-%v1)/(%t2-%t1)
                while %_1[%i]==. & %i<=_N {   /* goes back and fills */
                   qui replace %_1=%v1+%Slope*(_time-%t1) in %i
                   mac def i=%i+1
                }
             }
             mac def i=%i+1
         }
    macro shift
    }
}



* ===========  Substitution ===========================

* --------- parses INITIAL VALUES in the command option ----------
* if present, substitutes the model file values with the command values
capture confirm existence %_ival
if _rc==0 {
   parse "%_ival", parse(" ")
   while "%_1"!="" {
      qui replace %_1 in 1        /* assigns init. values, ex: ivpeso=10 */
      mac shift
   }
}

* ----------- parse PARAMETERS in the command option ---------
* if present, substitutes the model file values with the command values
capture confirm existence %_parv
if _rc==0 {
   parse "%_parv", parse(" ")
   while "%_1"!="" {
      capture qui replace  %_1         /* replace parameter value */
      if _rc!=0 { qui gen %_1 }
      mac shift
   }
}

* ----- Check for initial values -----
parse "%svar", parse(" ")
mac def flag=0
while "%_1"!="" {
   if %_1[1]==. {
        disp in red "Initial value not available for: %_1"
        mac def flag=1
   }
   macro shift
}
if %flag==1 {exit}


* ================= Computing ====================
disp in ye "Simulation:"
quietly {
mac def t=1                         /*  record number  */
while _time[%t]<%tend {
   mac def t0=%t-1
   disp in green "Time="_time[%t]

   * ----------- computes state equations ------------
   if %t>1 {
      parse "%svar", parse(" ")      /* list of state variables  */
      while "%_1"!="" {              /* %_1 holds, in turn, all state var. */
         mac def h2="seq%_1"
         replace bufs=%_1+%_dt*%%h2   in %t0
         replace %_1=bufs[%t0]        in %t
      mac shift
      }
   }
   * ----------- computes auxiliary equations --------
   parse "%avar", parse(" ")      /* list of aux. variables  */
   while "%_1"!="" {              /* %_1 holds, in turn, all aux var. */
        mac def h3="aeq%_1"
        replace %_1=%%h3  in %t
   mac shift
   }

   * ----------- computes rate equations -------------
   parse "%rvar", parse(" ")      /* list of rate variables  */
   while "%_1"!="" {              /* %_1 holds, in turn, all rate var. */
        mac def h0="req%_1"
        mac def code="fcod%_1"

        if %%code==0 {
           replace %_1=%%h0   in %t
        }

        else if %%code==1  {
           mac def h2 "idl%_1"
           mac def id=%t-%%h2[%t]/%_dt+0.5   /* information delay  */
           mac def id=int(%id)
           if %id>0 {
              replace bufr=%%h0       in %id
              replace %_1=bufr[%id]   in %t
           }
        }

        else if %%code==2 {
           mac def h1 "mdl%_1"
           mac def md=%t+%%h1[%t]/%_dt+0.5   /* material delay  */
           mac def md=int(%md)
           if %md<=_N {
              replace bufr=%%h0          in %t
              replace %_1=%_1+bufr[%t]   in %md
              * NOTE: the rate is added to a future rate, not substituted
           }
        }

        else if %%code==3 {
           mac def h1 "cvar%_1"
           mac def h2 "von%_1"
           mac def h3 "voff%_1"
           replace %_1=%%h0   in %t
           if %%h1[%t]<=%%h2[%t] | %%h1[%t]>=%%h3[%t] {replace %_1=0 in %t}
        }

        * if %id<1  the rate remains = 0
        * if %md>_N the rate is not calculated
   mac shift
   }

mac def t=%t+1

}
}
set more 0
/* %_options are all the graphic options specified by the user  */
gr %svar _time, %_options
drop bufs bufr zero
end