Revisiting Samuelson’s models, linear and nonlinear, stability conditions and oscillating dynamics

In this work, we reconsider the dynamics of a few versions of the classical Samuelson’s multiplier–accelerator model for national economy. First we recall that the classical one with constant governmental expenditure, represented by a linear second-order difference equation, is able to generate oscillations converging to the equilibrium for a wide range of values of the parameters, and give its analytic solution for all the possible cases. A delayed version proposed in the recent literature, represented by a linear third-order difference equation, is also considered. We show that also this model is able to produce converging oscillations, and give a complete analysis of the stability region of the equilibrium. A new simple nonlinear model is proposed, showing that it keeps oscillatory behavior, although coupled with other dynamics related to global effects. Our analysis confirms that the seminal work of Samuelson and simple modifications of it, may give powerful tools in the study of the business cycles.

We first recall that the original Samuelson's model, with constant governmental expenditure, which is described by a linear second-order difference equation, is capable to produce oscillations converging to the equilibrium value of the national income for different values of the parameters. We also consider the delayed version of Samuelson's model presented in Barros and Ortega (2019), which is described by a linear third-order difference equation. It also reveals its capability to produce oscillations converging to the equilibrium value, but for a smaller set of values for the parameters. Moreover, in order to describe a more realistic situation, we propose a simple nonlinear reformulation of the original Samuelson's model, with non constant governmental expenditure, and maintaining the delayed version presented in Barros and Ortega (2019). Our assumption on the governmental expenditure is that at each period it depends on some fixed costs plus a bounded quantity proportional to the income of the previous period. This simple nonlinear version of the Samuelson's model also evidences its ability to represent oscillations converging to the equilibrium value, although these may coexist with other global phenomena.
After this introduction, the rest of the paper is organized as follows. In Sect. 2 we reconsider the classical Samuelson's model and recall how its stability region is obtained depending on the two parameters, namely multiplier and accelerator, and give also its explicit analytic solution. In Sect. 3 we consider the delayed version of Samuelson's model presented in Barros and Ortega (2019). In that work the authors only determine the existence of the equilibrium point. Here we give the analytical description of its stability region, showing that it is in fact mainly associated with oscillatory dynamic behavior. In Sect. 3 we propose our simple nonlinear reformulation of the delayed version of Samuelson's model obtained introducing non constant governmental expenditure. The resulting nonlinear third-order difference equation has an attracting equilibrium for a wide region of values of the parameters, numerically detected. Some simulations show that the dynamics may converge to the equilibrium via oscillatory behavior, being accompanied by coexisting global phenomena which are related to the basin of attraction of the equilibrium. Section 4 concludes.

Original Samuelson's model
The original model by Samuelson is based on the following assumptions.
Assumption (A). National income T k at time k is given by the sum of three elements: consumption C k , private investment I k and governmental expenditure G k : (1) T k = C k + I k + G k Assumption (B). Consumption C k at time k depends on the past income and on marginal tendency to consume, modelled by a parameter c 1 , which is the multiplier coefficient, where 0 < c 1 < 1: Assumption (C). Private investment I k at time k depends on consumption changes and on the accelerator factor b, where b > 0 . Consequently, I k depends on national income changes: Assumption (D). Governmental expenditure G k at time k remains constant, G.
Hence, the national income is determined by the following second-order linear difference equation: which, by defining x k = T k , y k = T k+1 = x k+1 , leads to the 2D system whose fixed point is given by x = y = G 1−c 1 , that is: The stability is determined by the roots of the characteristic polynomial: which are given explicitly as follows: The stability occurs iff both eigenvalues are inside the unit circle. The stability region as a function of the two parameters is determined by the conditions P(1) = 1 − c 1 > 0, which is always satisfied, P(−1) = 1 + c 1 (1 + 2b) > 0, which is satisfied for b > − 1 2 − 1 2c 1 and P(0) < 1, which is satisfied for bc 1 < 1 . Thus, for any given 0 < c 1 < 1 and b > 0, the system is stable for 0 < b < 1 c 1 . Inside the stability region the eigenvalues may be real or complex. From the discriminant � = c 2 we have that the roots are complex conjugate when � < 0, that is, for any fixed value of c 1 ∈ (0, 1) the roots are complex conjugate for Notice that the two branches meet at c 1 = 1 giving b − (1) = b + (1) = 1, so that for any fixed value of c 1 ∈ (0, 1) at the bifurcation related to b = 1 c 1 the fixed point is a center (2) and for b > 1 c 1 it is a repelling focus. In Fig. 1 we show in color the stability region in the plane (c 1 , b), in yellow the region in which the fixed point is a stable focus, in orange the region in which it is a stable node.
Since the system in (4) is linear, we can write in explicit form its solution, in the considered range c 1 ∈ (0, 1) and b < 1 c 1 in which the fixed point is stable. As we know from elementary calculus of difference equations (see. e.g. Gumowsky andMira (1980), Elaydi (2005), only to cite a few) the second-order difference equation in (4) depends on two arbitrary constants which are determined by using the initial conditions, say T 0 and T 1 for k = 0 and k = 1 , respectively. Then, the solution is given for any k, i.e. T k for k ≥ 2, by the general solution of the homogeneous equation plus a particular solution, which in our case is simply the constant solution, i.e. the fixed point T * = G 1−c 1 . The solution of the homogeneous equation depends on the kind of eigenvalues, real or complex conjugate. When the eigenvalues s 1 and s 2 are real and distinct (i.e. for 0 < c 1 < 1 and 0 < b − (c 1 ) < 1 , orange region in Fig. 1), then the solution, for k ≥ 2, is given as follows: where Clearly, given the values of T k , for k ≥ 2, the other variables, which depend on T k , are easily obtained from their definitions: that is: In color the stability region of the fixed point T * = G 1−c1 is shown. In the yellow (resp. orange) region the eigenvalues are complex conjugate (resp. real) Tramontana and Gardini Economic Structures (2021) 10:9 from the initial conditions T 0 and T 1 for k = 0 and k = 1 , respectively, we have and When = 0 (i.e. for 0 < c 1 < 1 and b = b − (c 1 ) ) there are two coincident real roots and then the general solution T k , for k ≥ 2, becomes as follows: where Differently, when the eigenvalues s 1 and s 2 are complex conjugate ( � < 0 ), that is for having modulus √ c 1 b so that then the general solution T k , for k ≥ 2, which is spiraling around the fixed point, becomes as follows: where The other variables C k and I k , which depend on T k , are easily obtained from their definitions given in (12). This original system by Samuelson has been considered also recently in Ortega and Barros (2020), where the authors claim to give the explicit solutions for T k (and also C k and I k ). However, their solution is not correct, since it does not correspond to the one given above, and there is no discussion with respect to the eigenvalues, complex conjugate or real (distinct or coincident). We have so proved the following

Proposition 1
The Samuelson's model defined in (1), (2) and (3) with constant gov- The analytic solution of the equation given in (4) as a function of two arbitrary initial conditions T 0 and T 1 for k = 0 and k = 1, respectively, are given in (10), (11) when the eigenvalues (8) are real and distinct, in (16) and (17) when the eigenvalues are real and coincident, in (20) and (21) when the eigenvalues are complex conjugate. The other variables C k and I k are obtained from their definitions given in (12).

Delayed Samuelson's Model and stability analysis
The reformulated delayed version of Samuelson's model presented in Barros and Ortega (2019) is based on the following assumptions: Assumption (1) is equivalent to Assumption (A). National income T k at time k equals to the sum of three elements: consumption C k , private investment I k and governmental expenditures assumed constant G (as in Assumption (D)): Assumption (2). Consumption C k at time k is a linear function of the incomes of the two preceding periods: where c 1 and c 2 are positive constant, c 1 > 0 , c 2 > 0 , and 0 < c 1 + c 2 < 1 . The governmental expenditures G in the model by Barros and Ortega (2019) are included in the consumption C k , but since these are assumed constant, it does not differ from assuming them in the income T k .
Assumption (3) is equivalent to Assumption (C). Private investment I k at time k, depends on consumption changes and on the positive accelerator factor b > 0 : so that we get the following equation: Summarizing, by using (22), (23) and (25) the national income is determined via the following third-order linear difference equation: or, equivalently, (26) we have the fixed point, given by and under the assumption on the parameters (being 0 < c 1 + c 2 < 1) it is always positive and unique.
Notice that for c 1 + c 2 = 1 the fixed point becomes infinite and for c 1 + c 2 > 1 it is T * < 0, so that the condition c 1 + c 2 ≥ 1 can be considered unfeasible, and we assume the constraint 0 < c 1 + c 2 < 1, which, as we shall see, corresponds to the first condition for the stability region (associated with an eigenvalue equal to +1).
So it is now interesting to investigate the stability of the fixed point, as a function of the parameters (all positive), that is, under the given constraints: Let us introduce the three-dimensional map, by defining x k = T k , y k = T k+1 = x k+1 , z k = T k+2 = y k+1 , so that we have the 3D linear system and in the 3D space X = (x, y, z) the fixed point is The coefficient matrix is given by so that we recover the previous result, that is, from det (I − A) = 1 − (c 1 + c 2 ) we have the condition on the invertibility of (I − A) , given by c 1 + c 2 � = 1. The constant G is just a scaling factor, so it could be assumed G = 1 , and, as in the original Samuelson's model, it only influences the value of the fixed point T * .
The stability of the fixed point depends on the eigenvalues of the matrix A, that is, X * is attracting if all the eigenvalues (say (ξ 1 , ξ 2 , ξ 3 ) ) of A are smaller than 1 in absolute value. The eigenvalues are given by the roots of the characteristic polynomial: 33) a 1 = −c 1 (1 + b) = −(ξ 1 + ξ 2 + ξ 3 ) a 2 = bc 1 − c 2 (1 + b) = +(ξ 1 ξ 2 + ξ 1 ξ 3 + ξ 2 ξ 3 ) a 3 = bc 2 = −ξ 1 ξ 2 ξ 3 Tramontana and Gardini Economic Structures (2021) 10:9 In a recent work (Gardini et al. 2021), it is shown that assuming that P(0) = a 3 satisfies |a 3 | < 1, the stability conditions are as follows: In our case, the determinant is Det(A) = P(0) = a 3 = bc 2 > 0 so that it is always positive, and the condition reduces to bc 2 < 1. So we have to consider the two conditions 0 < c 2 < 1 b and 0 < c 2 < 1 − c 1 , that is: Substituting the expressions in (33) in conditions (34) we obtain: The equalities in the three conditions given in (36), coupled with (35), give the boundary of the stability region in the parameter space. Condition (i) is associated with an eigenvalue equal to 1, condition (ii) is associated with an eigenvalue equal to −1 , and condition (iii) with complex conjugate eigenvalues in modulus equal to 1. It follows that for any fixed value of the parameter b, b > 0 , we can have the boundaries of the stability region determined (besides c i > 0 ) by the conditions in (35) and where The boundaries of the region determined by the three conditions (i-iii) are given by curves in the parameter plane (c 1 , c 2 ) of equation: Let us first prove the following Proposition 2 For any b > 0 it holds that: Proof The straight lines c 2 = 1 b and (r 2 ) are intersecting in (c 1 , . The intersection point of the two straight lines (r 1 ) and (r 2 ) also comes soon by direct computation. Finally by direct computation we have 1 To comment the stability region of the fixed point, let us distinguish the two cases 0 < b ≤ 1 and b > 1, which have different properties in the parameter plane (c 1 , c 2 ).
The segment of (r 2 ) connecting the two points ( −1 1+2b , 0) and (0, 1 1+2b ) is inside the triangle bounded by the segment of c 2 = 1 − c 1 so that it is an active constraint. Notice that the point (0, 1 1+2b ) approaches (0, 1) as b tents to zero. (That is, increasing b from zero the stability region decreases, or, equivalently, decreasing b from 1 the stability region increases).
From the properties given in Proposition 2, the curve in the first quadrant of the parameter plane (c 1 , c 2 ) of equation c 2 = c 2 (c 1 ) is always external to the stability region, which is only bounded by segments of the straight lines (r 1 ) and (r 2 ). Some examples are shown in Fig. 2.
It follows that changing the parameters inside the stability region, either it is crossed the boundary at which the fixed point becomes infinite (segment of (r 1 ) ), or it is crossed the boundary at which one eigenvalue is equal to −1 (segment of (r 2 )).
(II) b > 1 Increasing b from 1, the curve of equation c 2 = c 2 moves in the opposite direction (since now it is 1 b < 1 ), intersecting the segment of (r 1 ) on the boundary of the stability region as long as it holds 1 < b < b * , examples are shown in Fig. 3a, b, while for b > b * the straight line (r 1 ) is external to the stability region, whose boundary is given by a segment of (r 2 ) and a segment of (r 3 ), an example is shown in Fig. 3c.
For the characteristic polynomial in (32) it is known (see e.g. Gardini et al. (2021)) that considering the discriminant then for D < 0 three real roots exist, at D = 0 two real roots merge and for D > 0 there are two complex conjugate roots and a real one. By using the expressions given in (33) (40) D = 4a 3 2 − a 2 1 a 2 2 + 4a 3 1 a 3 − 18a 1 a 2 a 3 + 27a 2 3 Tramontana and Gardini Economic Structures (2021) 10:9 we have numerically evaluated the discriminant D and shown the graph of D = 0 in Figs. 2, 3. As we can see, the stability region is mainly associated with a pair of complex eigenvalues, and for b > b * it is always so.
Notice that for a specific set of parameters we can numerically evaluate the eigenvalues and then, since the model is linear, we can easily write the explicit solution analytically as a functions of three initial conditions.
As an example, considering the values b = 1.5 , c 1 = 0.62, c 2 = 0.1 which is close to the boundary of (r 3 ), the fixed point T * = G 1−(c 1 +c 2 ) = 3.571429G is attracting, with one real negative eigenvalue and two complex conjugate ones, so that we have oscillations converging to the equilibrium, as shown in Fig. 4.
As a final remark, we can state that the delayed version of Samuelson's model reduces the stability region of the fixed point.

Nonlinear Samuelson's model
Let us come back to the original Samuelson's model, where the governmental expenditure may be not constant, G k , and combine this with the delayed version of Samuelson's  (c 1 , c 2 ). In a at b = 1 . In b at b = 0.8 . The black curve corresponds to the discriminant D = 0 given in (40). In the orange region a pair of complex eigenvalues exists. In the yellow region the three eigenvalues are real Fig. 3 In color the stability region of the fixed point T * for b > 1 is shown in the parameter plane (c 1 , c 2 ). In a at b = 1.2 , the black curve denotes the discriminant D = 0 given in (40). In b at b = b * . In c at b = 3 . In the orange region a pair of complex eigenvalues exist. In the yellow region the three eigenvalues are real model by Barros and Ortega (2019), by using G k in (22) in place of a constant value, keeping all the other assumptions as in the delayed model of the previous section.
The governmental expenditure G k may consist of a constant amount G (which includes, for example, fixed essential costs and other fixed expenses), and a part which depends on the income of the previous period, taking values in a bounded range, for example given by γ arctan ( T k−1 ), where γ is a parameter governing the width of the range, thus leading to our Assumption (4): This assumption is consistent, for instance, with the well known Wagner's hypothesis (Wagner 1890). According to this there is a causation that runs from economic growth to growth in government expenditure. This hypothesis has been tested empirically in several countries (see Gemmell 1990;Ansari 1993;Courakis et al. 1993;Dollery and Singh 1998;Peacock and Scott 2000, among the others). Thus, we get the following nonlinear equation leading to the three-dimensional nonlinear system The fixed point of the nonlinear system is now given by the solution of the following nonlinear equation: and it is easy to see that it can be considered as the intersection point of an increasing function ( γ arctan(T * ) ) with a straight line ( T * (1 − c 1 − c 2 ) − G ) with positive slope and negative offset, as qualitatively shown in Fig. 5 (where the function γ arctan(T * ) Thus, we may expect a unique fixed point in the positive side, as shown by the blue circle in Fig. 5, which is an intersection of the red curve with the blue straight line, but it is also possible to have three fixed points, as shown by the green circles in Fig. 5 which are the intersections of the red curve with the green straight line. However, this second situation may be considered unfeasible, since the two more fixed points must belong to the negative side.
Let us consider the positive fixed point. Since we cannot have its value analytically, the local stability analysis cannot be performed analytically. However, by using numerical simulations, we can see that the stability region is quite similar to the one of the linear case, as shown in Fig. 6a. 1 In that figure, at a value of b such that 0 < b < 1, we can see that the upper left boundary corresponds to a flip bifurcation (as in the linear case of the previous section). An example soon after the flip bifurcation is shown in Fig. 6b projected in the (x, y), where the unstable set of the fixed point converges to an attracting 2-cycle. In the upper right boundary the bifurcation is related to to the fixed point approaching infinity, similarly to the linear case of the delayed model considered in the previous section.
Also when b > 1 we have a stability region quite similar to the one observed in the linear case, as shown in Fig. 7a, with the upper left boundary leading to a flip bifurcation with a behavior similar to the one observed above (as shown in Fig. 7b), and the right side of the boundary of the stability region now includes also a border associated with a Neimark-Sacker bifurcation. It is worth noting that in the nonlinear model global bifurcations may occur, so that, also in the case of a unique fixed point, other attracting sets may appear. An example is shown in Fig. 8a, and enlarged in Fig. 8b where the attracting fixed point is close to the Neimark-Sacker bifurcation, but it coexists with another attractor, belonging to a closed invariant curve. From the enlargement we can argue that it is a cycle of very high period. The global effects involve regions in the phase space far from the equilibrium of interest, and it is relevant in order to determine the width of the basin of attraction of the equilibrium. The section of the basin of attraction on the plane z = T * is shown in Fig. 8c, and in the same figure we also have reported the projection of the two attractors. In it, yellow points converge to the fixed point while red points converge to the wide cycle, showing that the separator is probably a repelling closed invariant curve. As in fact, approaching the Neimark-Sacker bifurcation value, this curve shrinks to the fixed point in a subcritical bifurcation. Fig. 6 In a the yellow region represents the stability region of the fixed point T * in the parameter plane (c 1 , c 2 ) at b = 0.8 , G = 100, γ = 1 , initial condition (400, 400, 400). In b projection of the trajectory of the nonlinear model in the plane (x, y) = (T k , T k+1 ) showing the attracting 2-cycle existing at b = 0.8, c 1 = 0.2, c 2 = 0.584618 Fig. 7 In a the yellow region represents the stability region of the fixed point T * in the parameter plane (c 1 , c 2 ) at b = 1.5 , G = 100, γ = 1,initial condition (400, 400, 400). In b projection of the trajectory of the nonlinear model in the plane (x, y) = (T k , T k+1 ) showing the attracting 2-cycle existing at b = 1.5, c 1 = 0.2, c 2 = 0.45001

Conclusions
In the present work we have shown how strong is the ability of the Samuelson's model to represent stable oscillations converging to an equilibrium value. This has been shown in three cases. Namely, in Sect. 1 we have reconsidered the classical Samuelson's multiplier-accelerator model for national economy with constant governmental expenditure. The stability region has been analytically determined as well as its analytic solution in all the possible cases of real or complex eigenvalues. In Sect. 2 we have reconsidered the delayed version of Samuelson's model presented in Barros and Ortega (2019), performing the analysis of the stability region as a function of the three parameters, showing that indeed the stability region is mainly related to oscillatory behavior. In Sect. 3 we have proposed a simple nonlinear reformulation of the original Samuelson's model, with non constant governmental expenditure, maintaining the delayed version presented in Barros and Ortega (2019). Here we have performed a numerical investigation of the stability region, showing that it is quite similar to the one presented in Sect. 2, and also that the dynamics are mainly oscillating towards the stable equilibrium.
Our analysis shows that the seminal work of Samuelson (1939), and simple modifications of it, may give powerful tools in the study of the business cycles. A calibration of the parameters with realistic outcomes is left for future work.