Fortran vs VFI-Toolkit

Comparing Fortran and VFI-Toolkit for a finite-horizon life-cycle model with endogenous labor

Short summary: the VFI toolkit run times are not bad compared to Fortran. Value function iteration and distribution are about half an order of magnitude slower compared to Fortran. The calculation of age-dependent model statistics is the real bottleneck: too slow, especially considering that I set simoptions.whichstats to skip gini, lorenz and quintiles.

I have been working on a small replication exercise comparing a Fortran
implementation and a MATLAB/VFI-Toolkit implementation of the same
finite-horizon life-cycle model with endogenous labor supply.

The model has assets a, an exogenous Markov productivity shock \eta, age
j, and a permanent productivity type \theta. The household chooses
next-period assets a' and labor l. Labor is endogenous, but conditional on
a' it has a closed-form solution, so in the VFI-Toolkit implementation I do
not declare labor as a separate d decision grid. Instead the return function
uses the action-state inputs (a',a,\eta) and computes labor internally:

l(a',a,\eta,j,\theta) = \min\left\{\max\left\{ \nu + (1-\nu)\frac{a' - [(1+r)a+pen_j]}{w e_j \theta \eta}, 0\right\}, 1-10^{-10}\right\}

for working ages, and l=0 in retirement.

The replication code is here:

The Fortran code is in codes_fortran/. The MATLAB/VFI-Toolkit code is in
codes_matlab/. The model description and timing table are in tex_pdf/.

Runtime comparison

The table below reports the current runtime decomposition. The Fortran code is
compiled with Intel oneAPI ifx using /O2. The MATLAB code uses the
VFI-Toolkit finite-horizon routines with the toolkit’s GPU path when available.

Step Fortran (s) MATLAB/VFI-Toolkit (s) MATLAB/Fortran
solve_household / VFI 0.182225 0.896763 4.92x
Distribution 0.013436 0.077600 5.78x
Aggregation / model moments 0.001490 2.965766 1,990.45x
Total 0.197152 3.940129 19.99x

Overall, the toolkit implementation is about 20 times slower than the Fortran
implementation in this test. However, most of the total difference comes from
the calculation of model moments, not from the value-function iteration itself.

The VFI step is slower by about half an order of magnitude:

\log_{10}(0.896763/0.182225) \approx 0.69.

The distribution step is similar:

\log_{10}(0.077600/0.013436) \approx 0.76.

By contrast, the model-moment step is slower by more than three orders of
magnitude:

\log_{10}(2.965766/0.001490) \approx 3.30.

In levels, the aggregation/model-moment step accounts for about 79 percent of
the total MATLAB-minus-Fortran runtime gap in this run.

Issue with age-conditional standard deviations

One discrepancy I ran into concerns standard deviations conditional on age when
using permanent types. The toolkit output for grouped age profiles gave means
that looked fine, but the grouped standard deviations conditional on age did
not match the statistic I needed for the comparison with Fortran.

In particular, I wanted the standard deviation of an object x conditional on
age j, pooling across permanent types \theta_i. What worked was to use the
toolkit’s per-type age statistics and then combine them manually using the law
of total variance:

\operatorname{Var}(x \mid j) = \sum_i p_i \left[ \operatorname{Var}(x \mid j,\theta_i) + \left(E[x \mid j,\theta_i] - E[x \mid j]\right)^2 \right].

Then

\operatorname{SD}(x \mid j) = \sqrt{\operatorname{Var}(x \mid j)}.

This manual grouped-standard-deviation calculation is implemented in
codes_matlab/f_model_moments.m. It uses the toolkit-generated per-type means
and standard deviations, gathers the age statistics once, and then constructs
the grouped standard deviation outside the toolkit profile object.

I am posting this partly to document the comparison, and partly to ask whether
there is a better toolkit-native way to request this exact age-conditional,
permanent-type-pooled standard deviation.

3 Likes

This is a very nice example, thanks for sharing the code.

Is there a way to make the calculation of the moments faster and also to avoid the problem you discovered with aggregation of types?

2 Likes

Thanks @javier_fernan

I’m not sure if you can speed up the moment calculation in the toolkit: I already used whichstats to do only a subset of the available statistics. Of course you can reduce the “functions to evaluate” to be passed to the lifecycle moments function, but this means you don’t get the information you need.

My philosohy when using the toolkit: (this is a personal opinion, @robertdkirkby and others might have a different view).

Use the toolkit to do the policy functions and the distribution since that is the tricky part that takes most of the computing time. Moreover, programming it from scratch (in an efficient manner) takes time as well, even with AI. The toolkit is really optimized for this.
However, when it comes to compute, say, the coefficient of variation of some variable by age, do it yourself. It will be much faster than the toolkit and is relatively easy to program with AI. Quite often for a paper you want specific model moments calculated in a specific way and in the end is simpler to program this part youself. In this way you are also forced to check that the dimensions of the policy functions and the distribution are correct.

Later on, you can integrate your manual computations in the calibration routines of the toolkit by using fun_custom_stats.

Finally, if you look at my run times results, you can see that value function iteration and distribution are competitive with a hand-crafted Fortran code, while the age statistics are not.

2 Likes

First, a quick comment that the standard deviation bug with ptype models is now fixed:

1 Like