# Implementing a Stock-Flow Consistent macroeconomic model in Python

Over the last few years I’ve been developing, on and off, FSIC, a Python package for macroeconomic modelling. It’s a (perennial) work in progress. While I’m pretty happy with the user-facing elements, the implementation has always been a bit of a sandbox for me to learn Python. And it shows. It’s been a learning process and there’s a lot I’d do differently now.

This post documents a streamlined/barebones re-implementation of the core of the package with an example.

## Godley and Lavoie’s (2007) model of government money with portfolio choice

Chapter 4 of Godley and Lavoie (2007) sets out Model *PC*, an 11-equation
Stock-Flow Consistent (SFC)
model of government
money with portfolio choice. The model extends Model
*SIM* from Chapter 3 of
Godley and Lavoie (2007). Specifically, in addition to (high-powered) money,
Model *PC* introduces short-term government securities (bills), which are an
interest-paying financial asset. With two financial assets in the model,
households must make a portfolio decision to allocate their wealth between
money and bills, balancing the:

*transactions*motive for holding cash*speculative*motive for investing in financial assets which yield a return

The national income identity is the same as in Model *SIM*, reflecting a
closed, pure services economy. Production/income equals
the sum of consumption and government expenditure . The equation numbers match those in Godley and Lavoie (2007):

Bills held by households pay interest at the prevailing rate . Interest payments are an additional and taxable source of income for households and appear in the equations for both disposable income and taxes :

Each period, households receive interest on their holdings of bills from the previous period and at the previous period’s interest rate. Interest payments enter the above equations with a one-period lag.

To avoid complications to do with capital gains, the assumption is that bills have a fixed price of one unit.

# Portfolio allocation

In the model, households make a two-stage, hierarchical decision:

- a
*consumption*decision: how much to spend (and, conversely, how much to save) - a
*portfolio*decision: how to allocate savings (wealth)

Consumption is a function of disposable income and accumulated wealth . The and parameters are the propensities to consume out of each, respectively:

Wealth changes each period according to the difference between disposable income and consumption i.e. saving:

The portfolio decision involves a choice to allocate wealth between money and bills.

Godley and Lavoie (2007) follow the Brainard-Tobin approach in which households wish to hold a certain proportion of their wealth in bills and the rest as money. The proportions adjust for:

- the rate of return on bills: demand increases with the interest rate (the speculative motive)
- the ratio of disposable income to wealth (the transactions motive)

Tobin’s wealth/adding-up constraint holds such that the responsiveness to the interest rate of household demand for bills enters the money demand equation with the same parameter value but the opposite sign. The same is true for the parameter on the ratio of disposable income to wealth :

With only two financial assets, Equation 4.7 implies Equation 4.6A and vice versa. The model takes the demand for money as the residual after households decide on their demand for bills:

# Closing the model with endogenous money

The government deficit (change in debt, ) is the difference between outlays and revenues:

Outlays comprise expenditures and interest payments on outstanding debt. Revenues come from taxes as well as the profits the central bank returns to the government . These profits come from the interest earned on government debt held by the central bank. As the residual purchaser of bills, the central bank holds bills that households are not willing to hold at the given interest rate:

The central bank provides money to households on demand (money is endogenous) of an amount equal to the central bank’s additional holdings of bills:

The interest rate is exogenous as a consequence:

## Implementing the model in Python

The aim of FSIC is to allow the user to:

- specify a model in a convenient format which resembles the algebraic representation (as above)
- generate a Python class definition from the model specification
- use the newly-generated Python class to set values, solve the model and inspect the model results

# 1. Specification

The script for Model *PC* is below, as a Python string. Things to note:

- Scripts can include comments. As in the Python language, use
`#`

. - Time indexes follow the variable name in square brackets. Negative indexes
denote lags and positive indexes denote leads. For example,
`V[-1]`

refers to the one-period lag of wealth. - For contemporaneous variables, time indexes are optional:
`Y`

is the same as`Y[0]`

. - The current implementation does nothing special with them, but you can
denote parameters with braces e.g.
`{theta}`

. For equation errors, such as from an econometrically-estimated equation, you can use angle brackets e.g.`<epsilon>`

. - As in the Python language, use parentheses for implied line continuation, as
in the equation for
`Bh`

below.

# 2. Generation

The first step is to parse the script above to generate a list of
*symbols*. Symbols are the terms (variables and functions) that make up the
model. `Symbol`

objects are subclasses of the Python
`NamedTuple`

and store:

- the
`name`

of the symbol e.g.`Y`

- the
`type`

of the symbol; for example:`Y`

is an endogenous variable because there is an equation to determine its value`G`

is an exogenous variable because it only ever appears on the right-hand side of an equation`theta`

is a parameter

- the longest
`lags`

and`leads`

; for example:`Y`

has no lags or leads: only the contemporaneous values enter into the model each period`V`

has one lag from the wealth accumulation equation: a one-period lag in the model means that we can’t solve the very first period (because we’d always need a value for the period before)

- the
`equation`

representation (only for endogenous variables); a normalised representation with explicit time indexes - the
`code`

representation (again, only for endogenous variables), which is the Python code that enters the final model object

Parse the model script with
`parse_model()`

.

```
Symbol(name='Y', type=<Type.ENDOGENOUS: 3>, lags=0, leads=0, equation='Y[t] = C[t] + G[t]', code='self._Y[t] = self._C[t] + self._G[t]')
Symbol(name='G', type=<Type.EXOGENOUS: 2>, lags=0, leads=0, equation=None, code=None)
Symbol(name='theta', type=<Type.PARAMETER: 4>, lags=0, leads=0, equation=None, code=None)
Symbol(name='V', type=<Type.ENDOGENOUS: 3>, lags=-1, leads=0, equation='V[t] = V[t-1] + (YD[t] - C[t])', code='self._V[t] = self._V[t-1] + (self._YD[t] - self._C[t])')
```

My re-implementation tries to minimise third-party dependencies
(NumPy is the only one). However,
`pandas`

DataFrames are convenient to summarise the
model’s symbols.

name | type | lags | leads | equation | code | |
---|---|---|---|---|---|---|

0 | Y | 3 | 0 | 0 | Y[t] = C[t] + G[t] | self._Y[t] = self._C[t] + self._G[t] |

1 | C | 3 | 0 | 0 | C[t] = alpha_1[t] * YD[t] + alpha_2[t] * V[t-1] | self._C[t] = self._alpha_1[t] * self._YD[t] + self._alpha_2[t] * self._V[t-1] |

2 | G | 2 | 0 | 0 | None | None |

3 | YD | 3 | 0 | 0 | YD[t] = Y[t] - T[t] + r[t-1] * Bh[t-1] | self._YD[t] = self._Y[t] - self._T[t] + self._r[t-1] * self._Bh[t-1] |

4 | T | 3 | 0 | 0 | T[t] = theta[t] * (Y[t] + r[t-1] * Bh[t-1]) | self._T[t] = self._theta[t] * (self._Y[t] + self._r[t-1] * self._Bh[t-1]) |

5 | r | 3 | -1 | 0 | r[t] = r_bar[t] | self._r[t] = self._r_bar[t] |

6 | Bh | 3 | -1 | 0 | Bh[t] = V[t] * (lambda_0[t] + lambda_1[t] * r[t] - lambda_2[t] * (YD[t] / V[t])) | self._Bh[t] = self._V[t] * (self._lambda_0[t] + self._lambda_1[t] * self._r[t] - self._lambda_2[t] * (self._YD[t] / self._V[t])) |

7 | theta | 4 | 0 | 0 | None | None |

8 | alpha_1 | 4 | 0 | 0 | None | None |

9 | alpha_2 | 4 | 0 | 0 | None | None |

10 | V | 3 | -1 | 0 | V[t] = V[t-1] + (YD[t] - C[t]) | self._V[t] = self._V[t-1] + (self._YD[t] - self._C[t]) |

11 | lambda_0 | 4 | 0 | 0 | None | None |

12 | lambda_1 | 4 | 0 | 0 | None | None |

13 | lambda_2 | 4 | 0 | 0 | None | None |

14 | Hh | 3 | 0 | 0 | Hh[t] = V[t] - Bh[t] | self._Hh[t] = self._V[t] - self._Bh[t] |

15 | Bs | 3 | -1 | 0 | Bs[t] = Bs[t-1] + (G[t] + r[t-1] * Bs[t-1]) - (T[t] + r[t-1] * Bcb[t-1]) | self._Bs[t] = self._Bs[t-1] + (self._G[t] + self._r[t-1] * self._Bs[t-1]) - (self._T[t] + self._r[t-1] * self._Bcb[t-1]) |

16 | Bcb | 3 | -1 | 0 | Bcb[t] = Bs[t] - Bh[t] | self._Bcb[t] = self._Bs[t] - self._Bh[t] |

17 | Hs | 3 | -1 | 0 | Hs[t] = Hs[t-1] + Bcb[t] - Bcb[t-1] | self._Hs[t] = self._Hs[t-1] + self._Bcb[t] - self._Bcb[t-1] |

18 | r_bar | 2 | 0 | 0 | None | None |

Generate a Python class from the list of symbols with
`build_model()`

. This
returns a class with which to instantiate model objects. This class is derived
from an internal base class. The base class handles generic tasks like variable
access and model solution. The derived class adds model-specific information
such as the variable names and the code for the equations.

The automatically-generated underlying code for the derived class is in the
`CODE`

attribute.

```
class Model(BaseModel):
ENDOGENOUS = ['Y', 'C', 'YD', 'T', 'r', 'Bh', 'V', 'Hh', 'Bs', 'Bcb', 'Hs']
EXOGENOUS = ['G', 'r_bar']
PARAMETERS = ['theta', 'alpha_1', 'alpha_2', 'lambda_0', 'lambda_1', 'lambda_2']
ERRORS = []
LAGS = 1
LEADS = 0
def _evaluate(self, t):
self._Y[t] = self._C[t] + self._G[t]
self._C[t] = self._alpha_1[t] * self._YD[t] + self._alpha_2[t] * self._V[t-1]
self._YD[t] = self._Y[t] - self._T[t] + self._r[t-1] * self._Bh[t-1]
self._T[t] = self._theta[t] * (self._Y[t] + self._r[t-1] * self._Bh[t-1])
self._r[t] = self._r_bar[t]
self._Bh[t] = self._V[t] * (self._lambda_0[t] + self._lambda_1[t] * self._r[t] - self._lambda_2[t] * (self._YD[t] / self._V[t]))
self._V[t] = self._V[t-1] + (self._YD[t] - self._C[t])
self._Hh[t] = self._V[t] - self._Bh[t]
self._Bs[t] = self._Bs[t-1] + (self._G[t] + self._r[t-1] * self._Bs[t-1]) - (self._T[t] + self._r[t-1] * self._Bcb[t-1])
self._Bcb[t] = self._Bs[t] - self._Bh[t]
self._Hs[t] = self._Hs[t-1] + self._Bcb[t] - self._Bcb[t-1]
```

# 3. Solution

Initialise a new model instance by passing a sequence of time periods. This
represents the timespan of the model and could be, for example, a `range()`

or
a list of strings.

Under the bonnet, model variables are NumPy arrays. By default the arrays are
of type `float`

and have starting values of `0.0`

. Keyword arguments can
over-ride the default value e.g. `alpha_1=0.6`

. I do this below for the
behavioural parameters of the model i.e. the `alpha`

and `lambda`

variables.

After initialisation, you can get and set variable values by either:

- attribute e.g.
`pc.G = 20`

, which sets*all*values of`G`

to`20.0`

i.e. preserving the`float`

data type. Standard NumPy array indexing also works e.g.`pc.G[1:] = 20`

to set all but the first period’s value to`20.0`

- key e.g.
`pc['r_bar'] = 0.025`

or`pc['r_bar'][16:] = 0.035`

. If you pass a 2-tuple of indexes, the first index references the variable and the second index references a*named*period or range e.g.`pc['r_bar', 1960:] = 0.035`

sets to 3.5% from 1960 onwards. This saves having to keep track of array positions

Following Section 4.5.1 of Godley and Lavoie (2007) (‘The puzzling impact of interest rates’), the code below sets up the model to begin at its steady state at an interest rate of 2.5% and to simulate a step-change in the interest rate to 3.5% in 1960. Some values are from Gennaro Zezza’s EViews implementation of the same model.

As well as the model variables, other useful attributes of the model object include:

`span`

, which stores the periods of the model that were passed at initialisation`names`

, a list of strings of the model’s variable names`values`

, a two-dimensional NumPy array of the model variable values, with one row per variable and one column per period`status`

, a list of strings to record the solution state of each period as one of:`-`

: still to be solved (no attempt made)`.`

: solved successfully`F`

: failed to solve

`iterations`

, a list of integers to record the number of iterations to convergence each period (models are solved numerically); initialised to`-1`

Again, a `pandas`

DataFrame is a convenient way to inspect the model’s values.

Y | C | YD | T | r | Bh | V | Hh | Bs | Bcb | Hs | G | r_bar | theta | alpha_1 | alpha_2 | lambda_0 | lambda_1 | lambda_2 | status | iterations | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

1945 | 0.0 | 0.0 | 0.0 | 0.0 | 0.025 | 64.9 | 86.5 | 21.6 | 86.5 | 21.6 | 21.6 | 20.0 | 0.025 | 0.2 | 0.6 | 0.4 | 0.635 | 5.0 | 0.01 | - | -1 |

1946 | 0.0 | 0.0 | 0.0 | 0.0 | 0.000 | 64.9 | 86.5 | 21.6 | 86.5 | 21.6 | 21.6 | 20.0 | 0.025 | 0.2 | 0.6 | 0.4 | 0.635 | 5.0 | 0.01 | - | -1 |

1947 | 0.0 | 0.0 | 0.0 | 0.0 | 0.000 | 64.9 | 86.5 | 21.6 | 86.5 | 21.6 | 21.6 | 20.0 | 0.025 | 0.2 | 0.6 | 0.4 | 0.635 | 5.0 | 0.01 | - | -1 |

1948 | 0.0 | 0.0 | 0.0 | 0.0 | 0.000 | 64.9 | 86.5 | 21.6 | 86.5 | 21.6 | 21.6 | 20.0 | 0.025 | 0.2 | 0.6 | 0.4 | 0.635 | 5.0 | 0.01 | - | -1 |

1949 | 0.0 | 0.0 | 0.0 | 0.0 | 0.000 | 64.9 | 86.5 | 21.6 | 86.5 | 21.6 | 21.6 | 20.0 | 0.025 | 0.2 | 0.6 | 0.4 | 0.635 | 5.0 | 0.01 | - | -1 |

Solve the model using
`solve()`

. (See
the method’s docstring for information on various options.)

The full version of FSIC implements a graph algorithm to order the equations
for more efficient solution i.e. for convergence to occur in fewer
iterations. The barebones package I’m using here doesn’t include that algorithm
and I increase the maximum number of iterations (the `max_iter`

keyword
argument) to ensure convergence.

Y | C | YD | T | r | Bh | V | Hh | Bs | Bcb | Hs | G | r_bar | theta | alpha_1 | alpha_2 | lambda_0 | lambda_1 | lambda_2 | status | iterations | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

1945 | 0.0 | 0.0 | 0.0 | 0.0 | 0.025 | 64.9 | 86.5 | 21.6 | 86.5 | 21.6 | 21.6 | 20.0 | 0.025 | 0.2 | 0.6 | 0.4 | 0.635 | 5.0 | 0.01 | - | -1 |

1946 | 106.5 | 86.5 | 86.5 | 21.6 | 0.025 | 64.9 | 86.5 | 21.6 | 86.5 | 21.6 | 21.6 | 20.0 | 0.025 | 0.2 | 0.6 | 0.4 | 0.635 | 5.0 | 0.01 | . | 101 |

1947 | 106.5 | 86.5 | 86.5 | 21.6 | 0.025 | 64.9 | 86.5 | 21.6 | 86.5 | 21.6 | 21.6 | 20.0 | 0.025 | 0.2 | 0.6 | 0.4 | 0.635 | 5.0 | 0.01 | . | 101 |

1948 | 106.5 | 86.5 | 86.5 | 21.6 | 0.025 | 64.9 | 86.5 | 21.6 | 86.5 | 21.6 | 21.6 | 20.0 | 0.025 | 0.2 | 0.6 | 0.4 | 0.635 | 5.0 | 0.01 | . | 101 |

1949 | 106.5 | 86.5 | 86.5 | 21.6 | 0.025 | 64.9 | 86.5 | 21.6 | 86.5 | 21.6 | 21.6 | 20.0 | 0.025 | 0.2 | 0.6 | 0.4 | 0.635 | 5.0 | 0.01 | . | 101 |

Having solved the model, we can use matplotlib to recreate Figures 4.3 and 4.4 of Godley and Lavoie (2007).

The model shows how interest rates *increase* economic activity (the puzzling
impact) because there are no mechanisms for high interest rates to depress
aggregate demand.

See here for this post as a Jupyter notebook along with supporting Python code.

## References

Godley, W., Lavoie, M. (2007)
*Monetary economics: An integrated approach to
credit, money, income, production and wealth*,
Palgrave Macmillan

Zezza, G. (2006) ‘EViews macros for building models in *Wynne Godley and Marc
Lavoie* Monetary Economics: An Integrated Approach to Credit, Money, Income,
Production and Wealth’

http://gennaro.zezza.it/software/eviews/glch04.php

[Accessed 07/07/2018]