Part 2: Fixing the Stella R Python Script

Background

Welcome to my continued adventures in translating Stella models to R. To catch you up in case you haven’t been following along (i.e. a high probability that this applies to you) I’ve been retooling the StellaR.py script that I found on the internet. The latest version on the linked site is buggy and old, and written in Python 2. Additionally, I was given such a huge Stella model to work with that it would throw almost every exception in the book! I am almost finished with the new version of the StellaR.py conversion script (upgraded to Python 3), and though it’s not perfect (i.e. surely someone in the future will make some edits and thereby push the versioning forward), I’m happy to say it works: it translates old-version Stella models to R scripts (+1 for the open source movement). It does this by formatting the equation layer to be input into the ode() function in the package deSolve in R.

What I’m Currently Working On

In the old StellaR.py translation script, there was a note that the “delay function should be developed”…well, yes, it should! DELAY() was a function native to Stella, and the concept is really simple: it’s just a delayed version of whatever variable is being input. So, in Stella, let’s say the line is…

DELAY(X,1)

It’s just the variable X delayed by one time step, and the default initial value is 0. Practically in our model, it’s most commonly used to indicate e.g. the previous time step’s value alongside the current value in the system of equations. In an R data frame, this is a piece of cake, done with 1 line of indexing and a fill value at the beginning (I have a line just like it in several of my scripts).

The problem for me right now is conceptualizing how this is done in the context of the ode() function, which I am just at the beginner stage of learning. My big question about ode(), I guess, is “how does it work?” More specifically, is it iterating through time steps and calculating each value of each equation (and related, is that why it’s so slow? the code that runs it should be fast in terms of the languages (C/FORTRAN) they’re written in…)? Is this then generating a data frame as it goes as part of the object that results from the ode() function, or is it doing all of this “under the hood,” and if so, what hope do I have to index these variables as they’re being written? It feels like it should be an easy answer, and hopefully it will be, but it also seems like I’ll need to be able to access the values, and the resulting data frame, as it’s being calculated.

So I guess my biggest question is: can I reference prior values of the equations within/as the ode() is running, or do I need to handle it separately with the dede() function? My current line of thought is to write my own custom DELAY() function, which takes as arguments all of X’s specifications for the ode() function (e.g., the derivative equation, initial value, etc.) and just have it pass all of this to a dede() with the specified lag, returning the variable values at each time step. Here’s a pseudo-code jotting down of an example of my thoughts:

DELAY <- function(model_name,derivative,IV,lag,times) {
      lag_model <- <-function(t,Y,parms,...) {
      Time <<- t
      delayed <- model_name
      ddelayed <- derivative
      }
   dede(y=IV,times=time,func=lag_model)
}

T1_2previous_yeardry0 <- DELAY(T1_yeardry0,dT1_yeardry0,0,2,times)

A followup question though is whether or not this will run “alongside” the same time steps as the ode() function is taking, and thus return the values in step.

Leave a Reply

Your email address will not be published. Required fields are marked *