Appendix > Gekko solvers > Newton-Fair-Taylor

The Fair-Taylor (FT) algorithm can be perceived as performing Gauss-Seidel over time, instead of over equations. If the model is simulated over the period t1 to t2, this simulation is repeated again and again (where the values of any leaded variable are simply taken from the databank), until the leaded variable(s) converge (that is, do not move from iteration to iteration). The particular way each of the simulations from t1 to t2 are performed is irrelevant here, as long as the simulations solve the model.

Consider this one-equation model:

y = 0.1*y[-1] + 0.2*y + 0.3*y[+1] + 100

If the model is simulated from 2001 to 2004, the first simulation (2001) will take y[+1] as the 2002 value, the second simulation (2002) will take y[+1] as the 2003 value, the third simulation (2003) will take y[+1] as the 2004 value, and the last simulation (2004) will take y[+1] as the 2005 value. Regarding the 2005 value, there is the question of terminal conditions, but we might imagine that option 'exo' is used regarding terminals, so that the real (databank) value regarding y[2005] is used. We might consider this the first FT-iteration, that produced new values regarding y[2001], y[2002], y[2003] and y[2004].

In the next FT-iteration, consider the first sub-simulation, that is, the simulation of the year 2001. In this simulation, y[+1] is taken as the y[2002] value that was computed in the first FT-iteration. And likewise, regarding the simulation of the year 2002, y[+1] is taken as the y[2003] value that was computed in the first FT-iteration, and so on. The last simulation (2004) will still use the fixed value of y[2005], because we use the 'exo' option regarding terminals.

This process is repeated over and over, and in many cases it converges nicely. It can also be damped, and in spirit, the method is very similar to the Gauss-Seidel algorithm that solves the individual periods (hence, Fair-Taylor is also called Gauss-Seidel over time). Still, the process is by no means guaranteed to converge, and in cases with heavy intertemporal influences (for instance, if the parameter regarding y[+1] were 0.9 instead of 0.3), the 'signals' from the leaded variables from FT-iteration to FT-iteration may be slow to propagate from the last periods to first periods. If the intertemporal influence is too heavy, for instance with 0.9*y[+1] instead of 0.3*y[+1], the Fair-Taylor algorithm cannot solve the problem (no matter how much damping is used etc.).

In such cases, a more robust Newton solver is often used. The most common way to progress is typically to 'stack' the equations, eliminating all lags and leads. For instance, the system regarding 2001-2004 could be written in the following way:

y1 = 0.1*y0 + 0.2*y1 + 0.3*y2 + 100
y2 = 0.1*y1 + 0.2*y2 + 0.3*y3 + 100
y3 = 0.1*y2 + 0.2*y3 + 0.3*y4 + 100
y4 = 0.1*y3 + 0.2*y4 + 0.3*y5 + 100

This system is time-less (or static) in the sense that there are no lags or leads, but only variables with names that indicate the period. Provided that some values regarding y0 (= y[2000]) and y5 (= y[2005]) are provided, the system can be solved with a 'normal' solver (that is, a solver that solves one period at a time).

Now, such a system can become quite large, because the equations are unfolded (stacked) like in the example above. For a 100 year simulation period, the model becomes 100 times larger than a normal model, and the Jacobi matrix obtains 100*100 = 10.000 times more elements. There are a lot of tricks and methods to alleviate this explosion of variables, but implementing such tricks (for instance sparse matrices) is a time-consuming process.

Gekko implements such a stacked time solver, but it is mostly suited for smaller models, or for limited time periods. For larger models solved over long time periods, an alternative solver can be used. Provided that there are not too many variables with leads, this solver is quite powerful. The solver is called Newton-Fair-Taylor, and the basic principles of this solver is presented in the following section.

Newton-Fair-Taylor

This solver is activated by means of the following option:

option solve forward method = nfair; //default is 'fair'

The underlying idea is to view Fair-Taylor as an iterated process, where the goal is to find a fixed point in this process. If we limit the problem to one variable (y) containing leads, the process can be stated as follows:

ynew = F(yold)

So yold is a nx1 vector of initial values for y regarding the simulation period (with n observations), and this is fed to the normal solution algorithm, where leaded variables are just used at face value. The normal solution algorithm simulates the n periods, which produces n new values for y. The process is repeated over and over, and as soon as ynew = yold, the problem is solved.

This system can be linearized around some particular vector yold, by performing a small perturbation ε on one of the elements of yold, simulating the model over n periods, and observing the effects on ynew, compared to a simulation without ε. These perturbations are performed for each period, and produce n effects (for each of the periods contained in ynew). All in all, we obtain a n x n matrix of effects. For example, if the model y = 0.1*y[-1] + 0.2*y + 0.3*y[+1] is simulated over 4 periods, we obtain the following matrix:

1 2 3 4

1 0.0000 0.0000 0.0000 0.0000

2 0.3750 0.0469 0.0059 0.0007

3 0.0000 0.3750 0.0469 0.0059

4 0.0000 0.0000 0.3750 0.0469

The first row shows what happens to y[2001], y[2002], y[2003] and y[2004], if y[2001] is changed by one unit. In 2001, here is no effect of changing y[2001], since this only affects the starting values. As there are no effects in 2001, there are no effects in 2002-2004 either, so the first row is full of zeroes. The second row shows what happens, if y[2002] is changed by one unit. In 2001 this affects y, via the leaded variable. The effect is 0.3/(1-0.2) = 0.375. This in turn affects y in 2002 via the lag 0.1*y[-1]. This effect is 0.375*0.1/(1-0.2) = 0.0469, and so on.

This matrix is used to transform the error (ynew-yold) into a new guess regarding y. In order to do this, the 4x4 identity matrix is subtracted, and the resulting matrix is inverted. This matrix can then be combined with the error, to produce a new guess.

With option 'const' (option solve forward terminal = const) regarding terminal values, the idea is that regarding the terminal value y[2005], this value is set equal to y[2004]. Hence, the terminal value is assumed constant ('const') in relation to the y[2004] value, and in a stacked system, this would amount to the following:

y1 = 0.1*y0 + 0.2*y1 + 0.3*y2 + 100
y2 = 0.1*y1 + 0.2*y2 + 0.3*y3 + 100
y3 = 0.1*y2 + 0.2*y3 + 0.3*y4 + 100
y4 = 0.1*y3 + 0.2*y4 + 0.3*y4 + 100

Note the last equation, where y4 is used instead of y5. Using the 'const' option changes the incidence matrix, which would now be the following:

1 2 3 4

1 0.0000 0.0000 0.0000 0.0000

2 0.3750 0.0469 0.0059 0.0012

3 0.0000 0.3750 0.0469 0.0094

4 0.0000 0.0000 0.3750 0.0750

This changes the last column, for instance the element (4, 4) is changed from 0.0469 to 0.0750, that is a factor (1-0.2)/(1-0.2-0.3) = 1.6. With "OPTION solve forward terminal feed = internal", the incidence matrix will look like the matrix above, whereas with this option set to 'feed = external', the incidence matrix will look like the first one shown (with element (4, 4) = 0.0469).

With 'feed = internal', a linear model with leads will solve in one Fair-Taylor iteration, regardless of the coefficients (provided that the model is solvable). Using 'feed = external', more Fair-Taylor iterations are needed.

It would be tempting to provide a 'growth' option regarding terminal values, corresponding to this stacked system:

y1 = 0.1*y0 + 0.2*y1 + 0.3*y2 + 100
y2 = 0.1*y1 + 0.2*y2 + 0.3*y3 + 100
y3 = 0.1*y2 + 0.2*y3 + 0.3*y4 + 100
y4 = 0.1*y3 + 0.2*y4 + 0.3*y4*y4/y3 + 100

The last equation corresponds to y5/y4 = y4/y3, and the problem is that this equation may be hard to solve when solved on its own (because the right-hand side contains a term with y4 squared). Because of this, a 'growth' option regarding terminal values is not provided at the moment.

All in all, the Newton-Fair-Taylor solver ('nfair') is quite powerful regarding forward-looking models. If the number of variables with leads is limited, the incidence matrix does not become too large, and filling it up with coefficients does not become too time consuming. With k lead variables, the incidence matrix becomes (k*n) x (k*n) instead of n x n, and the matrix will contain cross-effects from one lead-variable to another. Such a matrix may be time-consuming to compute, since in principle k*n simulations are needed to fill it (not all simulations need to solve the full period, but there is still a lot of work to be done).

Looking at the original incidence matrix, an obvious idea springs to mind:

1 2 3 4

1 0.0000 0.0000 0.0000 0.0000

2 0.3750 0.0469 0.0059 0.0007

3 0.0000 0.3750 0.0469 0.0059

4 0.0000 0.0000 0.3750 0.0469

It is clear that the matrix is repetitive, with the same number (cf. the colors) running diagonally downwards. Using this idea, and assuming that the coefficients do not change too much over time in non-linear models, only k simulations are needed to fill the matrix (k being the number of leaded variables).

Some care must be taken regarding the coefficients corresponding to terminal values (if 'const' terminal condition is used), but a quite good approximation of the real matrix can probably be produced with little effort using this idea.

Newton-Fair-Taylor example

Provided the following model file:

------------- y.frm -----------------------------------

frml _i y = 0.1*y[-1] + 0.2*y + 0.3*y[+1] + 100;

-------------------------------------------------------

You may try the following statements:

RESET;
OPTION solve forward method = nfair;
OPTION solve forward dump = yes;
TIME 2001 2004;
MODEL y;
CREATE y;
SERIES <2000 2000> y = 200;
SIM;
PRT y;
PRT #ft_1;

Option 'nfair' is set regarding forward solving (default regarding terminals is always 'const'). Option 'dump' is set, so that the gradient matrix can be shown afterwards (#ft_1). The result is the following:

+++ NOTE: There are 1 variable(s) with leads: Fair-Taylor algorithm is used

#1: Gauss simulation 2001-2004 took 0.00 sec -- 10/14/8.8 iterations (min/max/avg)

Gradient 1 of 4 (var 1 per 2001)

Gradient 2 of 4 (var 1 per 2002)

Gradient 3 of 4 (var 1 per 2003)

Gradient 4 of 4 (var 1 per 2004)

#2: Gauss simulation 2001-2004 took 0.01 sec -- 10/14/43.4 iterations (min/max/avg)

Newton-Fair-Taylor (leads) algorithm converged in 3 NFT-iterations (1.34 sec)

y %

2001 243.4254 21.71

2002 249.1343 2.35

2003 249.8830 0.30

2004 249.9766 0.04

#ft_1

1 2 3 4

1 0.0000 0.0000 0.0000 0.0000

2 0.3750 0.0469 0.0059 0.0012

3 0.0000 0.3750 0.0469 0.0094

4 0.0000 0.0000 0.3750 0.0750

The gradient matrix uses four simulations (one for each period), and the algorithm converges in 3 Newton-Fair-Taylor iterations, which is always the case for a linear model. Note that in this simulation, the actual databank value of y in 2005 is not used at any point (it is missing value anyway). From the results, it is seen that y is developing towards it's equilibrium value 250. This value can be found by removing all lags/leads in the equation y = 0.1*y[-1] + 0.2*y + 0.3*y[+1] + 100, so that it becomes y = 0.1*y + 0.2*y + 0.3*y + 100 (which has y = 250 as solution).

Setting "OPTION solve forward terminal = exo", and setting y[2005] = 200 would change the solution completely:

2001 242.2783 21.14

2002 246.0754 1.57

2003 242.1082 -1.61

2004 230.2635 -4.89

In this case, the y[2004] is 'attracted' to y[2005] = 200, and hence cannot attain a value close to 250 anymore. So the choice of terminal condition should not be taken too lightly, and the 'const' option makes much more sense than the 'exo' option in most cases.