# Clear the environment
rm(list=ls(all=TRUE))
# Set number of periods
Q=200
# Set number of scenarios (including baseline)
S=4
# Set period in which shock/shift will occur
s=50
# Create (S x Q)-matrices that will contain the simulated data
Y=matrix(data=1,nrow=S,ncol=Q) # Income/output
C=matrix(data=1,nrow=S,ncol=Q) # Consumption
I=matrix(data=1,nrow=S,ncol=Q) # Investment
# Set fixed parameter values
c0=0
c1=0.75
I0=10
# Set baseline parameter values
Yc=matrix(data=80,nrow=S,ncol=Q) # output ceiling
D=matrix(data=10,nrow=S,ncol=Q) # investment floor (depreciation)
beta=matrix(data=1.25,nrow=S,ncol=Q) # accelerator coefficient
# Set parameter values for different scenarios
Yc[2,s:Q]=100 # scenario 2: increase in output ceiling
D[3,s:Q]=1 # scenario 3: increase in investment floor
beta[4,s:Q]=0.8 # scenario 4: reduction in accelerator coefficient
# Simulate the model by looping over Q time periods for S different scenarios
for (i in 1:S){
for (t in 2:Q){
for (iterations in 1:500){ # run the model 500-times in each period
#Model equations
#(1) Output
Y[i,t] = min(C[i,t]+I[i,t], Yc[i,t])
#(2) Consumption
C[i,t] = c0 + c1*Y[i,t-1]
#(3) Investment
I[i,t] = I0 + max(beta[i,t]*(Y[i,t-1] - Y[i,t-2]), -D[i,t])
} # close iterations loop
} # close time loop
} # close scenario loop
16 A Hicksian Business Cycle Model with Floors and Ceilings
Overview
This model captures some key features of John R. Hicks (1950)’ nonlinear model of business cycle fluctuations. The model builds on Samuelson (1939)’s linear multiplier-accelerator model analysed in Chapter 2 and Chapter 9, but gives it a different twist. While Samuelson (1939)’s business cycle model requires repeated external shocks to produce sustained cycles, Hicks (1950)’s model produces permanent, shock-independent fluctuations that are bounded by a floor and a ceiling. The goods market equilibrium is assumed to be locally unstable due to a strong accelerator effect, i.e. investment reacts strongly to increases in aggregate output, which creates more output and boosts investment even further. However, aggregate output is constrained by a ceiling that may stem from full employment (and a fixed labour supply). On the downward trajectory, output is constrained by a floor that is given by a minimum amount of investment needed to replace the depreciated capital stock.
We present a version of the model that is adapted from Gabisch and Lorenz (1989), chapter 2.1.1
The Model
\[ Y_{t}=\min[C_t + I_t, Y^C] \tag{16.1}\]
\[ C_{t} = c_0 + c_1 Y_{t-1}, \quad c_1 \in (0,1) \tag{16.2}\]
\[ I_{t} = I_0 + \max[\beta(Y_{t-1} - Y_{t-2}), -D], \quad \beta > 1, D >0 \tag{16.3}\]
where \(Y_t\), \(C_t\), \(I_t\), \(Y^C\), and \(D\), represent output, consumption, investment, the (full employment) output ceiling, and depreciation, respectively.
Equation 16.1 specifies aggregate income. For values below the ceiling \(Y^C\), aggregate income is determined by aggregate demand (consisting only of consumption and net investment). As soon as aggregate demand hits the ceiling, aggregate income is given by the ceiling. Equation 16.2 is a simple linear Keynesian consumption function. Equation 16.3 is the investment function. For investment levels above depreciation \(D\), net investment is given by an autonomous component \(I_0\) and an accelerator component that is responsive to the change in aggregate income.2 The reaction coefficient \(\beta\) is assumed to be sufficiently large to generate explosive dynamics. Once the change in aggregate income has fallen below the negative of the level of depreciation \(D\), net investment is given by the floor \(I_0 - D\).
Simulation
Parameterisation
Table 1 reports the parameterisation used in the simulation. Besides a baseline (scenario 1), we consider a fall in the output ceiling \(Y^C\) in scenario 2, an increase in the investment floor \(I_0 - D\) in scenario 3, and finally a decrease in the accelerator coefficient \(\beta\) (scenario 4).
Table 1: Parameterisation
Scenario | \(c_0\) | \(c_1\) | \(I_0\) | \(\beta\) | \(Y^C\) | \(D\) |
---|---|---|---|---|---|---|
1: baseline | 0 | 0.75 | 10 | 1.25 | 80 | 10 |
2: fall in output ceiling (\(Y^C\)) | 0 | 0.75 | 10 | 1.25 | 100 | 10 |
3: increase in investment floor (\(I_0 - D\)) | 0 | 0.75 | 10 | 1.25 | 80 | 1 |
4: decrease in accelerator coefficient (\(\beta\)) | 0 | 0.75 | 10 | 0.8 | 80 | 10 |
Simulation code
import numpy as np
# Set number of periods
= 200
Q
# Set number of scenarios (including baseline)
= 4
S
# Set period in which shock/shift will occur
= 50
s
# Create (S x Q)-arrays that will contain the simulated data
= np.ones((S, Q)) # Income/output
Y = np.ones((S, Q)) # Consumption
C = np.ones((S, Q)) # Investment
I
# Set fixed parameter values
= 0
c0 = 0.75
c1 = 10
I0
# Set baseline parameter values
= np.full((S, Q), 80) # output ceiling
Yc = np.full((S, Q), 10) # investment floor (depreciation)
D = np.full((S, Q), 1.25) # accelerator coefficient
beta
# Set parameter values for different scenarios
1, s:Q] = 100 # scenario 2: increase in output ceiling
Yc[2, s:Q] = 1 # scenario 3: increase in investment floor
D[3, s:Q] = 0.8 # scenario 4: reduction in accelerator coefficient
beta[
# Simulate the model by looping over Q time periods for S different scenarios
for i in range(S):
for t in range(1, Q):
for iterations in range(500): # run the model 500-times in each period
# Model equations
# (1) Output
= min(C[i, t] + I[i, t], Yc[i, t])
Y[i, t]
# (2) Consumption
= c0 + c1 * Y[i, t - 1]
C[i, t]
# (3) Investment
= I0 + max(beta[i, t] * (Y[i, t - 1] - Y[i, t - 2]), -D[i, t])
I[i, t]
# Close iterations loop
# Close time loop
# Close scenario loop
Plots
Figure 16.1 displays the model’s dynamics under the baseline scenario. The model generates cycles in aggregate output, consumption and investment that are constrained by the output ceiling \(Y^C=80\) and the investment floor \(I_0-D=0\). These cycles are permanent, i.e. they don’t require any external shocks, and they are (roughly) periodic, with a cycle length from peak (trough) to peak (trough) of around 13.5 periods. What generates the turning points? Suppose income is on an increasing trajectory. The accelerator effect on investment will amplify this process as firms will increase their investment, which raises aggregate demand and increases income further. When the output ceiling is reached, the change in income becomes zero, and investment demand decreases relative to the boom phase. This decrease in investment lowers aggregate output, and the accelerator effect again amplifies this tendency that now plays out on a downward trajectory until net investment hits the lower bound. Once investment has stopped its downward trajectory, aggregate income stops falling as well, which allows the accelerator effect to turn positive again and thereby induce a new cycle.
# Set start and end periods for plots
Tmax=100
Tmin =10
#Income, consumption, and investment
plot(Y[1, Tmin:(Tmax)],type="l", col=1, lwd=2, lty=1, xlab="", ylab="Y, C, I", ylim=range(I[1, Tmin:Tmax],Y[1, Tmin:(Tmax)]))
title(main="Income, consumption, and investment (baseline)", xlab = 'Time',cex=0.8,line=2)
lines(C[1, Tmin:Tmax],lty=2)
lines(I[1, Tmin:Tmax],lty=3)
legend("topleft", legend=c("Y", "C", "I"),
lty=1:3, cex=0.8, bty = "n", y.intersp=0.8)
In Figure 16.2, the output ceiling rises to \(Y^C = 100\) in \(t=50\). This raises the amplitude of the cycles. Interestingly, the ceiling is no longer binding, but the cycles still do not explode. This illustrates a feature of the model that has been noted before: it actually does not need a ceiling to generate bounded cycles (see, e.g., Puu, Gardini, and Sushko (2005)). The acceleration of income during a boom will become smaller and smaller, and eventually turn negative, thereby lowering investment and inducing a switch from boom to bust. Hence, cycles (and thus turning points) will occur even without floor and ceilings, but to prevent these cycles to become larger and larger, a floor is sufficient.
# Scenario 2: increase in ceiling
plot(Y[2, Tmin:(Tmax)],type="l", col=1, lwd=2, lty=1, xlab="", ylab="Y, C, I", ylim=range(I[2, Tmin:Tmax],Y[2, Tmin:(Tmax)]))
title(main="Scenario 2: Increase in output ceiling", xlab = 'Time',cex=0.8,line=2)
lines(C[2, Tmin:Tmax],lty=2)
lines(I[2, Tmin:Tmax],lty=3)
legend("topright", legend=c("Y", "C", "I"),
lty=1:3, cex=0.8, bty = "n", y.intersp=0.8)
An increase in the investment floor to \(I_0 - D = 9\) in period \(t=50\), displayed in Figure 16.3, has the effect of reducing the amplitude of the fluctuations. Intuitively, the milder the bust, the milder the recovery.
# Scenario 3: increase in investment floor
plot(Y[3, Tmin:(Tmax)],type="l", col=1, lwd=2, lty=1, xlab="", ylab="Y, C, I", ylim=range(I[3, Tmin:Tmax],Y[3, Tmin:(Tmax)]))
title(main="Scenario 3: Increase in investment floor", xlab = 'Time',cex=0.8,line=2)
lines(C[3, Tmin:Tmax],lty=2)
lines(I[3, Tmin:Tmax],lty=3)
legend("topright", legend=c("Y", "C", "I"),
lty=1:3, cex=0.8, bty = "n", y.intersp=0.8)
Finally, Figure 16.4 displays the scenario where the accelerator coefficient \(\beta\) falls, rendering investment much less sensitive to changes in aggregate income. As a result, the model’s regular cycles turn into damped oscillations that eventually converge to a stable equilibrium.
# Scenario 4: Reduction in accelerator coefficient
plot(Y[4, Tmin:(Tmax)],type="l", col=1, lwd=2, lty=1, xlab="", ylab="Y, C, I", ylim=range(I[4, Tmin:Tmax],Y[4, Tmin:(Tmax)]))
title(main="Scenario 4: Reduction in accelerator coefficient", xlab = 'Time',cex=0.8,line=2)
lines(C[4, Tmin:Tmax],lty=2)
lines(I[4, Tmin:Tmax],lty=3)
legend("topright", legend=c("Y", "C", "I"),
lty=1:3, cex=0.8, bty = "n", y.intersp=0.8)
### Plots (here only for baseline)
import matplotlib.pyplot as plt
# Set start and end periods for plots
= 100
Tmax = 10
Tmin
# Plot income (Y), consumption (C), and investment (I)
range(Tmin, Tmax), Y[0, Tmin:Tmax], label='Y', color='black', linewidth=2, linestyle='-')
plt.plot("Income, consumption, and investment (baseline)", fontsize=10)
plt.title('Time')
plt.xlabel('Y, C, I')
plt.ylabel(range(Tmin, Tmax), C[0, Tmin:Tmax], label='C', linestyle='--', color='black',)
plt.plot(range(Tmin, Tmax), I[0, Tmin:Tmax], label='I', linestyle=':', color='black',)
plt.plot(
# Add a legend to the top-left corner
='upper left', fontsize=8, frameon=False)
plt.legend(loc
# Show the plot
plt.show()
Directed graph
Another perspective on the model’s properties is provided by its directed graph. A directed graph consists of a set of nodes that represent the variables of the model. Nodes are connected by directed edges. An edge directed from a node \(x_1\) to node \(x_2\) indicates a causal impact of \(x_1\) on \(x_2\).
## Create directed graph
# Construct auxiliary Jacobian matrix for 5 variables:
# endogenous: (1) Y, (2) C, (3) I
# exogenous: (4) I0-D, (5) Yc
#Y C I I0-D Yc
M_mat=matrix(c(0,1,1, 0, 1, # Y
1,0,0, 0, 0, # C
1,0,0, 1, 0, # I
0,0,0, 0, 0, # I0-D
0,0,0, 0, 0),# Yc
5, 5, byrow=TRUE)
# Create adjacency matrix from transpose of auxiliary Jacobian and add column names
A_mat=t(M_mat)
# Create directed graph from adjacency matrix
library(igraph)
dg=graph_from_adjacency_matrix(A_mat, mode="directed", weighted= NULL)
# Define node labels
V(dg)$name=c("Y", "C", "I", expression(I[0]-D), expression(Y^C))
# Plot directed graph matrix
plot(dg, main="Directed graph of Hicks model", vertex.size=40, vertex.color="lightblue",
vertex.label.color="black", edge.arrow.size=0.3, edge.width=1.1, edge.size=1.2,
edge.arrow.width=1.2, edge.color="black", vertex.label.cex=1.2,
vertex.frame.color="NA", margin=-0.08)
### Directed graph
import networkx as nx
# Construct auxiliary Jacobian matrix for 5 variables:
# endogenous: (1) Y, (2) C, (3) I
# exogenous: (4) I0-D, (5) Yc
# Create the matrix M_mat
= np.array([[0, 1, 1, 0, 1], # Y
M_mat 1, 0, 0, 0, 0], # C
[1, 0, 0, 1, 0], # I
[0, 0, 0, 0, 0], # I0-D
[0, 0, 0, 0, 0]]) # Yc
[
# Create adjacency matrix from transpose of auxiliary Jacobian
= M_mat.T
A_mat
# Create directed graph from adjacency matrix using networkx
= nx.DiGraph(A_mat)
G
# Define node labels
= {0: "Y", 1: "C", 2: "I", 3: "$I_0 - D$", 4: "$Y^C$"}
nodelabs
# Plot directed graph
= nx.spring_layout(G, k=0.08)
pos =200, node_color="lightblue",
nx.draw_networkx(G, pos, node_size="black", width=1.2, arrowsize=10,
edge_color='->', font_size=8, font_color="black",
arrowstyle=True, labels=nodelabs)
with_labels"off")
plt.axis("Directed Graph of Lewis Model")
plt.title( plt.show()
It can thus be seen that the key exogenous variables in the Hicks (1950) model are the output ceiling \(Y^C\) and the investment floor \(I_0 - D\), which constrain the endogenous cycle generated by the interaction between the endogenous variables income \(Y\), investment \(I\), and consumption \(C\).
Analytical discussion
The dynamic properties of the Hicks (1950) model when the floor and ceiling are not binding are similar, but not entirely equivalent to the Samuelson (1939) model analysed in Chapter 9.
First, define \(X_t = Y_{t-1}\) to re-write the system as:
\[ X_{t}=\min[C_{t-1} + I_{t-1}, Y^C] \tag{16.4}\]
\[ C_{t} = c_0 + c_1 X_{t} \tag{16.5}\]
\[ I_{t} = I_0 + \max[\beta(X_{t} - X_{t-1}), -D]. \tag{16.6}\]
Then shift Equation 16.5 one period back and substite into Equation 16.4 to obtain:
\[ X_{t}=\min [c_0 + c_1 X_{t-1} + I_{t-1}, Y^C]. \tag{16.7}\]
Next, substitute Equation 16.7 into Equation 16.6:
\[ I_{t} = I_0 + \max[\beta( c_0 + (c_1 -1) X_{t-1} + I_{t-1}), -D]. \tag{16.8}\]
Equation 16.7 and Equation 16.8 are a two-dimensional dynamic system in \(X_t\) and \(It\).
For \(Y_t < Y^C\) and \(I_t > I_0 - D\), the system’s Jacobian matrix is given by:
\[ J=\begin{bmatrix} c_1 & 1 \\ \beta(c_1-1) & \beta \end{bmatrix}. \]
The characteristic polynomial yielding the eigenvalues of the Jacobian is:
\[\lambda^2-\lambda(c_1+\beta)+\beta=0,\]
where \(c_1+\beta =tr(J)\) and \(\beta = det(J)\).
Thus, we have:
\[ \lambda_{1,2} = \frac{(c_1+\beta) \pm \sqrt{(c_1+\beta)^2-4\beta}}{2}, \]
and two eigenvalues will be a pair of complex conjugates if \((c_1+\beta)^2-4\beta <0\) or \(\left( \frac{c_1+\beta}{2} \right)^2 < \beta\).
To assess stability, we need to check whether the modulus of the complex eigenvalue \(R\) is inside the unit circle, i.e. \(R \in (-1,1)\). For a two-dimensional system with complex eigenvalues, we have \(R=\sqrt{det(J)}\). Thus, the stability condition for the Hicks model is \(\beta <1\).
The following code generates a plot that displays the condition for cycles, \(\left( \frac{c_1+\beta}{2} \right)^2 < \beta\), and the stability condition, \(\beta <1\), in the \((\beta, c_1)\)-space:3
# Create function for cycle condition using beta as argument
cyc= function (beta) {
2*sqrt(beta)-beta
}
# Plot the two conditions in (beta, c1)-space
curve(cyc, from = 0, to = 4.5, col = 1, xlab=expression(beta), ylab=expression(c[1]) , main="",
lwd=1.5, n=10000, ylim=range(0,1.5))
abline(v=1, col=2)
legend("topright", legend = c("cycle condition", "stability condition"),
col = c(1, 2), lwd = 2)
# Create function for cycle condition using beta as an argument
def cyc(beta):
return 2 * np.sqrt(beta) - beta
# Define the range for beta
= np.linspace(0, 4.5, 10000)
beta_values
# Compute the corresponding c1 values using the cyc function
= cyc(beta_values)
c1_values
# Plot the cycle condition in (beta, c1)-space
='black', linewidth=1.5, label='cycle condition')
plt.plot(beta_values, c1_values, color
# Plot the stability condition as a vertical line at beta = 1
=1, color='red', linewidth=2, label='stability condition')
plt.axvline(x
# Set labels and title
r'$\beta$', fontsize=12)
plt.xlabel(r'$c_1$', fontsize=12)
plt.ylabel('')
plt.title(
# Set y-axis limit
0, 1.5)
plt.ylim(
# Add legend
='upper right')
plt.legend(loc
# Show the plot
plt.show()
The Hicks model requires combinations of \(c_1\) and \(\beta\) below the cycle condition and to the right of the stability condition to yield cycles that would be explosive without the floor.
We can also compute the eigenvalues and check the general stability conditions numerically.
### Stability analysis
## Baseline
# Construct Jacobian matrix
J_base=matrix(c(c1, 1,
beta[1,1]*(c1-1), beta[1,1]),
2, 2, byrow=TRUE)
# Obtain eigenvalues
ev_base=eigen(J_base)
(evals_base = ev_base$values)
[1] 1+0.5i 1-0.5i
# Obtain determinant and trace
(det=det(J_base)) # determinant
[1] 1.25
[1] 1.118034
#Check general stability conditions
print(1+tr+det>0)
[1] TRUE
print(1-tr+det>0)
[1] TRUE
print(1-det>0)
[1] FALSE
# Check analytical stability condition
beta[1,1] < 1
[1] FALSE
## Compare with scenario 4
J_scen4=matrix(c(c1, 1,
beta[4,s]*(c1-1), beta[4,s]),
2, 2, byrow=TRUE)
# Obtain eigenvalues
ev_scen4=eigen(J_scen4)
(evals_scen4 = ev_scen4$values)
[1] 0.775+0.4465143i 0.775-0.4465143i
# Obtain determinant and trace
(det=det(J_scen4)) # determinant
[1] 0.8
[1] 0.8944272
#Check general stability conditions
print(1+tr+det>0)
[1] TRUE
print(1-tr+det>0)
[1] TRUE
print(1-det>0)
[1] TRUE
# Check analytical stability condition
beta[4,s] < 1
[1] TRUE
For the baseline, the modulus of the complex eigenvalue is indeed larger than unity and the system thus unstable. This is due to the accelerator coefficient \(\beta\) being larger than unity. By contrast, in scenario 4, we have \(\beta < 1\) and the system is thus stable.
### Check cycle condition and compute cycle length (baseline)
# Check analytical cycle condition
(c1+beta[1,1])^2 - 4*beta[1,1] < 0
[1] TRUE
# Save real and imaginary part of complex eigenvalue
re=Re(evals_base[1])
im=Im(evals_base[1])
# Calculate cycle length
L=(2*pi)/(acos(re/mod_base))
L
[1] 13.55164
The analytical cycle condition confirms that the eigenvalues of the system will be complex and thus generate cycles. The implied cycle length is around 13.5 periods.
### Stability Analysis - Baseline
# Construct the Jacobian matrix for the baseline scenario
= np.array([[c1, 1],
J_base 0, 0] * (c1 - 1), beta[0, 0]]])
[beta[
# Obtain eigenvalues
= np.linalg.eig(J_base)
evals_base, _ print("Eigenvalues (Baseline):", evals_base)
# Obtain determinant and trace
= np.linalg.det(J_base) # determinant
det_base = np.trace(J_base) # trace
tr_base
# Calculate modulus of the first eigenvalue
= np.abs(evals_base[0])
mod_base print("Modulus (Baseline):", mod_base)
# Check general stability conditions
print("Stability conditions (Baseline):")
print(1 + tr_base + det_base > 0)
print(1 - tr_base + det_base > 0)
print(1 - det_base > 0)
# Check analytical stability condition
print("Analytical stability condition (Baseline):", beta[0, 0] < 1)
### Compare with Scenario 4
# Construct the Jacobian matrix for scenario 4
= np.array([[c1, 1],
J_scen4 3, s] * (c1 - 1), beta[3, s]]])
[beta[
# Obtain eigenvalues
= np.linalg.eig(J_scen4)
evals_scen4, _
# Print eigenvalues
print("Eigenvalues (Scenario 4):", evals_scen4)
# Obtain determinant and trace
= np.linalg.det(J_scen4) # determinant
det_scen4 = np.trace(J_scen4) # trace
tr_scen4
# Calculate modulus of the first eigenvalue
= np.abs(evals_scen4[0])
mod_scen4 print("Modulus (Scenario 4):", mod_scen4)
# Check general stability conditions
print("Stability conditions (Scenario 4):")
print(1 + tr_scen4 + det_scen4 > 0)
print(1 - tr_scen4 + det_scen4 > 0)
print(1 - det_scen4 > 0)
# Check analytical stability condition
print("Analytical stability condition (Scenario 4):", beta[3, s] < 1)
### Check cycle condition (baseline)
= (c1 + beta[0, 0])**2 - 4 * beta[0, 0] < 0
cycle_condition print("Cycle condition (Baseline):", cycle_condition)
# Save real and imaginary parts of the complex eigenvalue
= np.real(evals_base[0])
re = np.imag(evals_base[0])
im
# Calculate cycle length (L)
= (2 * np.pi) / np.arccos(re / mod_base)
L print("Cycle length (Baseline):", L)
References
Following Hommes (1995), we consider a stationary version of the model that abstracts from the growth path assumed in Hicks (1950).↩︎
This differs from Samuelson (1939), who assumed that investment reacts to changes in consumption, see Chapter 2.↩︎
To create the function for the cycle condition, we set \(\left( \frac{c_1+\beta}{2} \right)^2 = \beta\) and solve for \(c_1\), where we only consider positive values for \(c_1\) and \(\beta\).↩︎