Notes on using F2PY to call Fortran code from Python.

The application is the power series approximation of the Leontief inverse matrix, from input-output analysis in economics. The explanation and example below draw heavily on Miller and Blair (2009) and Almon (2017).

Application

Notation and accounting

In input-output analysis we define:

  • an matrix of inter-industry purchases (intermediate consumption):
  • an vector of primary inputs (value added):
  • an vector of final demand:
  • an vector of gross output:

By the accounting that underpins an input-output table:

where is a summation vector: an vector of ones.

The Leontief inverse matrix

Under the assumption of a Leontief production function, which exhibits constant returns to scale, we define , an matrix of technical (or input-output) coefficients with elements:

where is what industry (product) requires of industry (product) to produce one unit of output.

We can then express gross output, , as:

Solving for :

where is the (here, Type I) Leontief inverse matrix:

and is the identity matrix.

Equation 5 provides a way to answer the question, ‘How much output would be required from across the economy to support a certain amount of final demand?’. For a hypothetical vector of final demand, :

Power series approximation

Historically, matrix inversion was computationally expensive and, even today, may not be a good idea. As an alternative, we can think of the Leontief inverse as the sum of a power series:

where , etc.

With this, we can approximate a solution to Equation 5 as:

for some large enough for to be not significantly different from zero. This is the power series approximation we want to implement. Miller and Blair (2009) provide more detail on the derivation and, as well as Almon (2017), the proof by which this series converges for an input-output matrix, .

The power series also has a convenient economic interpretation. Each term in the series represents a round-by-round output effect to support a given amount of final demand. To supply our dinner, the butcher and the brewer must buy from the farmer; the baker must buy from the miller, who must also buy from the farmer etc.

Implementation

F2PY provides various ways to wrap Fortran code. Here, I do it the quick way but with a Python function to wrap the call to the compiled Fortran code. An abridged version of the Fortran subroutine is below:

subroutine power_series_approximation(f, A, q, m)
  implicit none

  integer, intent(in) :: m
  real(8), dimension(m), intent(in) :: f
  real(8), dimension(m, m), intent(in) :: A
  real(8), dimension(m), intent(out) :: q
  real(8), dimension(m) :: term
  integer :: i

  ! Copy `f` to initialise the series
  q = f
  term = f

  do i = 1, 25
     ! Calculate the next term in the series and add it to `q`
     term = matmul(A, term)
     q = q + term

     ! Test for convergence
     if(sum(term * term) < 0.000001) then
        exit
     end if

  end do

end subroutine power_series_approximation

Niceties the code above could do with:

  • control over the maximum number of iterations, rather than a hardcoded limit of 25
  • control over the convergence threshold, rather than a hardcoded value of 0.000001
  • some indication of whether the procedure converged

The full code implements the above here.

We then compile the code from the command line with f2py (here, the code above is in a file with name ‘leontief.f95’):

$ f2py -c leontief.f95 -m leontief

The Python wrapper (again, abridged) is as follows:

def power_method_fortran(f, A):
    from leontief import power_series_approximation

    q = power_series_approximation(f, A)
    return q.reshape(f.shape)

We can now call power_method_fortran() just like any other Python function.

The full code uses the example of the Tiny model from Almon (2017). The ‘power_series_example.py’ script compares the speed of three approaches to calculate from and :

  1. matrix inversion and then multiplication by Equations 5 and 6 above, in Python
  2. the power method, implemented in Python
  3. the power method, implemented in Fortran (as above)

All rely on NumPy to some degree.

Running the script on my machine, I get the following screen output:

$ python power_series_example.py
Times (speed gain relative to Python power method):
 - Python (matrix inversion): 2.38s (2.91x)
 - Python (power method): 6.92s (1.00x)
 - Fortran (power method): 0.28s (24.30x)

This shows a reasonable speed gain from doing this. In a one-off piece of analysis, the time savings probably aren’t worth it, though. You’re unlikely to make the coding (and debugging) time back. However, if you have an application that involves repeatedly calculating gross output (say, in a large-scale macroeconometric input-output model with changing input-output coefficients), it could be worth it.

See here for the full code and example.

References

Almon, C. (2017) The craft of economic modeling, Third, enlarged edition, Inforum
http://www.inforum.umd.edu/papers/TheCraft.html

Miller, R. E., Blair, P. D. (2009) Input-output analysis: Foundations and extensions, Second edition, Cambridge University Press