Monetary Economics in R: Model SIM

I've always been interested in macroeconomics so I started reading Godley and Lavoie's Monetary Economics book, it describes economies with increasingly complex models and I thought I could have fun implementing them in R.

In this post we will explore the first, and simplest, model: model SIM.

Model SIM

Model SIM is made of only eleven equations, I won't explain what they mean because I assume you are familiar with the book. (And if you're not, you should read it! Highly recommended).

\begin{aligned}
C_s & = C_d\\
G_s & = G_d\\
T_s & = T_d\\
N_s & = N_d\\
YD & = W \cdot N_s - T_s\\
T_d & = \theta \cdot W \cdot N_s\\
C_d & = \alpha_1 \cdot YD + \alpha_2 \cdot H_{h-1}\\
\Delta H_s & = H_s - H_{s-1} = G_d - T_d\\
\Delta H_h & = H_h - H_{h-1} = YD - C_d\\
Y & = C_s + G_s\\
N_d & = \frac{Y}{W}
\end{aligned}

where \theta < 1 and 0 < \alpha_2 < \alpha_1 < 1.

The SIM function

We can now implement the model in R. I decided to write it as a function that takes the exogenous variables as parameters and returns a data frame containing the simulation data.

#' Generate the data for the SIM model.
#' 
#' @param G governement expenditures
#' @param alpha consumption parameters
#' @param theta personal income tax rate
#' @param W nominal wage rate
#' @param N the number of iterations
#' @param initial.state the state at the "beginnig of the world"
#' @return the simulation result
SIM <- function(G = 20, alpha = c(0.6, 0.4), theta = 0.2, W = 1,
                N = 200, initial.state = 0) {
  # Calculate steady state values once
  Y.s <- G / theta
  YD.s <- (G * (1 - theta)) / theta
  # Initialize the data frame
  data <- data.frame(Year = rep(0, N),
                     C    = rep(0, N),
                     G    = rep(0, N),
                     t    = rep(0, N),
                     N    = rep(0, N),
                     YD   = rep(0, N),
                     DHs  = rep(0, N),
                     DHh  = rep(0, N),
                     Y    = rep(0, N),
                     H    = rep(0, N),
                     V.t  = rep(0, N),
                     Y.s  = rep(Y.s, N),
                     YD.s = rep(YD.s, N))
  data[1, ] = initial.state
  start <- data[1,]$Year
  # Update the data frame with the model values for each iteration
  for (i in 2:N) {
    Year <- start + i
    H_1 <- data[i-1,]$H
    Y <- (G + alpha[2] * H_1) / (1 - alpha[1] * (1 - theta))
    N <- Y / W
    t <- theta * Y
    YD <- Y - t
    DHs <- G - t
    C <- alpha[1] * YD + alpha[2] * H_1
    DHh <- YD - C
    H <- DHs + H_1
    V.t <- YD * (1 - alpha[1]) / alpha[2]
    data[i,1:11] <- c(Year, C, G, t, N, YD, DHs, DHh, Y, H, V.t)
  }
  return(data)
}

We can test it using the parameters provided in the book and check the result against Table 3.4.

res <- SIM()
head(res, 3)
#   Year        C  G        t        N       YD      DHs      DHh        Y        H Y.s
# 1    0  0.00000  0 0.000000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000   0
# 2    2 18.46154 20 7.692308 38.46154 30.76923 12.30769 12.30769 38.46154 12.30769 100
# 3    3 27.92899 20 9.585799 47.92899 38.34320 10.41420 10.41420 47.92899 22.72189 100
tail(res, 1)
#     Year  C  G  t   N YD          DHs          DHh   Y  H Y.s
# 200  200 80 20 20 100 80 5.329071e-14 5.684342e-14 100 80 100

Plotting model SIM

The data we generate is right, but viewing it as a table is not really informative or interesting.

Using the ggplot2 library we can generate 4 plots that looks like the ones in the book.

First we start by generating Figure 3.2 and Figure 3.3:

data <- SIM() # We use the default parameters
# Figure 3.2
fig.3.2 <- ggplot(data, aes(x=Year)) +
  labs(x="", y="") +
  theme(legend.direction="horizontal", legend.position="top") +
  scale_color_hue("") +
  xlim(0, 50) +
  geom_line(aes(y=YD, colour="Disposable Income YD")) +
  geom_line(aes(y=C, colour="Consumption C"))
print(fig.3.2)
# Figure 3.3
fig.3.3 <- ggplot(data, aes(x=Year)) +
  labs(x="", y="") +
  theme(legend.direction="horizontal", legend.position="top") +
  scale_color_hue("") +
  xlim(0, 50) +
  geom_line(aes(y=H, color="Wealth Level")) +
  geom_line(aes(y=DHh, color="Household Saving"))
print(fig.3.3)

This code will produce the two following images

Figure 3.2 Figure 3.3

Figure 3.1 and Figure 3.4 require us to increase the government spending by an amount ∆G. To do this we need the data for a stationary state and then shock it by increasing the governement expenditures G.

data.with.gov.expenses.increase <- function() {
  d <- SIM(G = 20)
  d <- tail(d, 10)
  d <- rbind(d, SIM(G = 25, initial.state = d[10,]))
  return(d)
}
data <- data.with.gov.expenses.increase()
# Figure 3.1
fig.3.1 <- ggplot(data, aes(x=Year)) +
  labs(x="", y="") +
  theme(legend.direction="horizontal", legend.position="top") +
  scale_color_hue("") +
  xlim(190, 250) +
  geom_line(aes(y=Y, colour="Income Y")) +
  geom_line(aes(y=Y.s, colour="Steady-State solution Y*"))
print(fig.3.1)
# Figure 3.4
fig.3.4 <- ggplot(data, aes(x=Year)) +
  labs(x="", y="") +
  theme(legend.direction="horizontal", legend.position="top") +
  scale_color_hue("") +
  xlim(190, 250) +
  geom_line(aes(y=V.t, color="Target Wealth")) +
  geom_line(aes(y=YD, color="Disposable Income")) +
  geom_line(aes(y=C, color="Consumption")) +
  geom_line(aes(y=H, color="Wealth"))
print(fig.3.4)

The figures generated by this snippet are shown below

Figure 3.1 Figure 3.4

In the last picuter we can't see the Target Wealth line because in this specific case, that is when \frac{1-\alpha_1}{\alpha_2} = 1, it's equal to Disposable Income.