Integration of economic model into existing framework

[twallema]

Currently, I have a working implementation of an economic production model in the notebook twallema-uncoupled-economic-model.ipynb on my branch twallema/uncoupled_economic_model. I would like to incorporate this model as much as possible using the existing modeling framework of Stijn and Joris. The economic model is in essence, a discrete-time system which can be written as,

def integrate(t, old_states, parameters)
    1. old_states --> new_states
return new_states

My economic model has 7 states, which are stratified in 63 economical sectors,

• x : productivity of sector i
• c : household consumption of good i
• f : exogeneous consumption of good i
• d: total demand for good i
• l : labor income in sector i
• O : total orders made by sector i
• S : stock matrix, per sector i, how much stock of good from sector j remain

All states are 63-dimensional vectors, except the last state, S, which is a 63x63 matrix. I have so far tried to use a manual definition of the state sizes but I find this complicates the code and is not a modular solution of the problem. Because reshape is used, this will only work if some states are square matrices of size [stratification_size * stratification_size]. The state S does not necessarily need to be returned to the end-user. It is not an interesting state but it is needed to calculate the new values of the other states. Maybe the issue can be circumvented by somehow initializing an initial stock matrix S in a place/way so that integrate can load the previous value, perform its calculation of the other states, and then store the new stock matrix S someplace (and repeat the above in an iterative manner).

Something like this,

def integrate(t, old_states, parameters)
    1. load old S from someplace
    2. old_states --> new_states
    3. store new S in someplace
return new_states

However, I need some pointers or inspiration on how to accomplish this.
@jorisvandenbossche @stijnvanhoey

[stijnvanhoey]

I just wanted to check the concept, which requires alterations on 2 levels:

• getting the initial conditions fit the framework (1)
• do a different reshape for the 1D versus the 2D variables to link scipy solver with out integrate function (2)
• do a different reshape for the 1D versus the 2D variables of the output (3)

A dummy model setup, just to verify the concept of handling the variable dimensions:

state_2d = "y2"    # y2 is a 2D variable, y1 is a 1D variable
state_names = ["y1", "y2"]
stratification_size = 3    # stratification level of 3

split_point = (len(state_names) - 1) * stratification_size

# User would write the integrate function and make sure it works for the envisioned (2D) dimensions
def integrate(t, y1, y2):
    """y1 1x3 ; y2 3x3  """
    dy1 = -0.3*y1
    dy2 = 0.1*y2
    return (dy1, dy2)

# `create_fun` method of the base class should make sure the function returned should do the following
def func(t, y, pars={}):
    """As used by scipy -> flattend in, flattend out"""
    
    # incoming y -> different reshape for 1D vs 2D variables  (2)
    y_1d, y_2d = np.split(y, [split_point])
    y_1d = y_1d.reshape(((len(state_names) - 1), stratification_size)),
    y_2d = y_2d.reshape((stratification_size, stratification_size))
    
    dstates = integrate(t, *y_1d, y_2d)
    return np.concatenate([np.array(state).flatten() for state in dstates])

# Setup/run model
initial_states = {"y1": [42., 42., 42.], 
                         "y2": np.ones((3, 3))}
time = [0, 5]
initial_states[state_2d] = initial_states[state_2d].flatten()  # Handle intial condition (1)
result = solve_ivp(func, time, list(itertools.chain(*initial_states.values())))

# To output: different reshape for 1D vs 2D variables to collect output  (3)
output = result["y"]
y_1d, y_2d = np.split(output, [split_point])
y_1d, y_2d.reshape((len(result["t"]), stratification_size, stratification_size))

Note: cfr. remarks of #178 about the usage of the scipy solver. Most of the reshaping (bookkeeping dimensions) is to fit it into the scipy solver. If in practice, the discrete version is always used (and the called method is an update of the states as such), this kind of bookkeeping is probably not required as the solve_discrete can pass whatever we want/define it to do.