Macro simulations
NOTE: The present notebook is coded in R. It relies heavily on the tidyverse ecosystem of packages. We load the tidyverse below as a prerequisite for the rest of the notebook - along with a few other libraries.
\(\rightarrow\) Don’t forget that code flows sequentially. A random chunk may not work if the previous have have not been executed.
library(tidyverse) # Package for data wrangling
library(readxl) # Package to import MS Excel files
library(latex2exp) # Package for LaTeX expressions
library(quantmod) # Package for stock data extraction
library(highcharter) # Package for reactive plots
library(patchwork) # Package for plot layout
library(plotly) # Library for interactive plots (3D!)The content of the notebook is based on the Macro Simulation website and the neoclassical case in particular.
1 Context
Until now, we have seen a lot of theory, but limited simulations… This notebook is intended on providing some elements in this direction. And to show if the models achieve convergence - or not.
2 The model
2.1 Production & co.
We first start with a Cobb-Douglas production function: \[Y=AK^aN^{1-a}, \quad a \in(0,1), \tag{1}\] where \(N\) is labour/employment/active population. Next, the real wages \[w=(1-a)A(K/N)^a \tag{2}\] This comes from the labour demand of profit-maximizing firms. Indeed, the profit is \[\Pi(N)=AK^aN^{1-a}-wN-\text{(other costs)}, \tag{3}\]
hence the FOC is \(\frac{\partial \Pi}{\partial N}=(1-a)A(K/N)^a-w=0\) - QED.
Government expenditures are tied by: \(G=T+B\), i.e., taxes plus debt. But the future value is \(G^f=T^f-(1+r)B\) because the interest of the debt must be reimbursed (no default here!). Solving for \(B\) gives
\[G+\frac{G^f}{1+r}=T+\frac{T^f}{1+r}. \tag{4}\]
Labor supply is given by \[N = 1-b_1/w. \tag{5}\]
One way to come to this result is the following and stems from arbitrage with respect to leisure. We can assume that the agent in the economy has utility \[U(C,N)= C+b \log(1-N),\]
where \(1-N\) represents (properly normalized) the time dedicated to leisure. The parameter \(b\) tunes how much leisure matters for the agent. At the same time, the agent needs to work to consume: \(C=wN\). There are thus two competing effects: increase work to consume more or decrease work to spend more time on leisure. The logarithmic function is key to obtaining the desired form.
The FOC with respect to N (replacing \(C\) in the above expression) is thus \[\frac{\partial U}{\partial N}=w-\frac{b}{1-N}=0 \quad \Leftrightarrow \quad w(1-N)=b \quad \Leftrightarrow \quad N=1-b/w\]
2.2 Consumption
In fact the true utility is more complex: \[U(C)=\log(C+b_1\log(1-N))+b_2\log(C^f),\] where \(C^f\) is future consumption (expected and treated as exogenous), which is subject to a budget constraint involving future output. The current constraint is
\[C=Y-T-S,\] i.e., what is consumed is what is produced, minus taxes minus what is saved at rate \(r\). The future version is \[C^f=Y^f-T^f+(1+r)S.\] Solving for \(S\),
\[C^f=Y^f-T^f+(1+r)(Y-T-C),\] where we can replace taxes by government spending thanks to Equation 4: \[C^f=Y^f-G^f+(1+r)(Y-G-C).\]
Hence, the utility becomes
\[U(C,N)=\log(C+b_1\log(1-N))+b_2\log(Y^f-G^f+(1+r)(Y-G-C))\]
Then,
\[\frac{\partial U}{\partial C}=\frac{1}{C+b_1\log(1-N)}-\frac{b_2(1+r)}{Y^f-G^f+(1+r)(Y-G-C)},\] so the FOC with respect to \(C\) is
\[C^f=b_2(1+r)(C+b_1\log(1-N)) \quad \Leftrightarrow \quad C=\frac{C^f}{b_2(1+r)}-b_1\log(1-N) \tag{6}\]
Also, the derivative with respect to \(N\) requires to recall \(Y= \text{some terms} + wN\):
\[\frac{\partial U}{\partial N}=\frac{-b_1}{(1-N)(C+b_1\log(1-N))}+\frac{b_2(1+r)w}{C^f}\] where we recognize a similar term from Equation 6, which gives
\[C=w\frac{(1-N)(C+b_1\log(1-N))}{b_1}-b_1\log(1-N)\]
and rearranging \[C+b_1\log(1-N)=w\frac{(1-N)(C+b_1\log(1-N))}{b_1} \Leftrightarrow b_1/w = (1-N)\] QED (same form as before!).
3 The code
In fact, the original code has some unnecessary parts if we are only interested in production or consumption.
Below, we keep only the required rows.
3.1 Initializing variables
S <- 5 # Number of scenarios (including baseline)
nb_sim <- 20 # Number of simulations
# Create and parameterise exogenous variables/parameters that will be shifted
G0 <- vector(length=S) # government expenditures
A <- vector(length=S) # productivity
Yf <- vector(length=S) # expected future income
b1 <- vector(length=S) # household preference for leisure
G0[] <- 1
A[] <- 2
Yf[] <- 1
b1[] <- 0.4# Set parameter values for different scenarios
G0[2] <- 2 # scenario 3: fiscal expansion
A[3] <- 3.5 # scenario 4: productivity boost
Yf[4] <- 0.2 # scenario 5: lower expected future income
b1[5] <- 0.8 # scenario 6: increased preference for leisure
#Set constant parameter values
a <- 0.3 # Capital elasticity of output
b2 <- 0.9 # discount rate
b3 <- 0.6 # household preference for money
K <- 5 # Exogenous capital stock
pe <- 0.02 # Expected rate of inflation
Gf <- 1 # Future government spending
# Initialise endogenous variables at arbitrary positive value
w <- matrix(1, ncol = S, nrow = nb_sim)
C = I = Y = r = N = P = w3.2 Running the loop
for (i in 1:S){
for (t in 2:nb_sim){
# Model equations
Y[t, i] <- A[i]*(K^a)*N[t-1, i]^(1-a) #(1) Cobb-Douglas production function
w[t, i] <- A[i]*(1-a)*(K^a)*N[t-1, i]^(-a) #(2) Labor demand
N[t, i] <- 1 - (b1[i])/w[t-1, i] #(3) Labor supply
#(4) Consumption demand
C[t, i] <- (1/(1+b2+b3))*(Y[t-1, i] - G0[i] + (Yf[i]-Gf)/(1+r[t-1, i]) -b1[i]*(b2+b3)*log(b1[i]/w[t-1, i]))
r[t, i] <- (I[t-1, i]^(a-1))*a*A[i]*N[t,i]^(1-a) #(5) Investment demand, solved for r
I[t, i] <- Y[t, i] - C[t,i] - G0[i] #(6) Goods market equilibrium condition, solved for I
}
}Let’s plot the result!
bind_cols(t = 1:nb_sim, as_tibble(Y)) |>
pivot_longer(-t, names_to = "scenario", values_to = "production") |>
ggplot(aes(x = t, y = production, color = scenario)) + geom_line() +
theme_classic() +
theme(legend.position = "top",
axis.title.x = element_blank())
Some scenarios are redundant (initial choices have no impact).
Let’s have a look at consumption…
bind_cols(t = 1:nb_sim, as_tibble(C)) |>
pivot_longer(-t, names_to = "scenario", values_to = "consumption") |>
ggplot(aes(x = t, y = consumption, color = scenario)) + geom_line() +
theme_classic() +
theme(legend.position = "top",
axis.title.x = element_blank())
Ok, some curves stop early… what happened? At time \(t=3\), the variable \(I\) became negative because \(G\) is too large!
What if we add a little noise on the production side? With small shocks to labor supply.
for (i in 1:S){
for (t in 2:nb_sim){
# Model equations
Y[t, i] <- A[i]*(K^a)*N[t-1, i]^(1-a) # (1) Cobb-Douglas production function
w[t, i] <- A[i]*(1-a)*(K^a)*N[t-1, i]^(-a) # (2) Labor demand
N[t, i] <- 1 - (b1[i])/w[t-1, i] + rnorm(1, sd = 0.02) # (3) Labor supply
}
}
bind_cols(t = 1:nb_sim, as_tibble(Y)) |>
pivot_longer(-t, names_to = "scenario", values_to = "production") |>
ggplot(aes(x = t, y = production, color = scenario)) + geom_line() +
theme_classic() +
theme(legend.position = "top",
axis.title.x = element_blank())
Ok, so we still have convergence (mostly).
What if the shocks are ten times larger?
set.seed(42)
for (i in 1:S){
for (t in 2:nb_sim){
# Model equations
Y[t, i] <- A[i]*(K^a)*N[t-1, i]^(1-a) # (1) Cobb-Douglas production function
w[t, i] <- A[i]*(1-a)*(K^a)*N[t-1, i]^(-a) # (2) Labor demand
N[t, i] <- 1 - (b1[i])/w[t-1, i] + rnorm(1, sd = 0.2) # (3) Labor supply
}
}
bind_cols(t = 1:nb_sim, as_tibble(Y)) |>
pivot_longer(-t, names_to = "scenario", values_to = "production") |>
ggplot(aes(x = t, y = production, color = scenario)) + geom_line() +
theme_classic() +
theme(legend.position = "top",
axis.title.x = element_blank())
Rocky, but still relatively consistent.
4 Calibration
4.1 Introduction
Modern economic models are complex. The dominating paradigm is that of dynamic stochastic general equilibrium (DSGE) models. The problem is that their complexity raises issues with respect to parameter estimation. If we follow Estimating DSGE Models: Recent Advances and Future Challenges, then we write \[\textbf{Y}_t=f(\textbf{S}_t,\textbf{V}_t, \boldsymbol{\theta}),\]
where
- \(\textbf{Y}_t\) is a vector of observables (GDP, inflation, etc.);
- \(\textbf{S}_t\) is a vector of state variables (we could assume they are observable for simplicity);
- \(\textbf{V}_t\) is a vector of shocks;
- \(\boldsymbol{\theta}\) is a vector of parameters.
This makes a lot! The aim is to find reasonable (plausible) values for \(\theta\). As discussed in the aforementioned paper, there are many ways to do so. But let’s here take a simple route.
4.2 Toy example
Let’s recycle the model we played with above, but with no shocks.
There is no scenario, parameters belong to the function, as described earlier.
model <- function(A, a, b1, K, nb_sim){
Y <- 1
w <- 0.5
N <- 0.5
for (t in 2:nb_sim){
Y[t] <- A * (K^a) * N[t-1]^(1-a) # (1) Cobb-Douglas production function
w[t] <- A * (1-a) * (K^a) * N[t-1]^(-a) # (2) Labor demand
N[t] <- 1 - b1 / w[t-1] # (3) Labor supply
}
cbind(Y, w, N)#[nb_sim]
}
model(1, 0.5, 0.3, 1, 15) Y w N
[1,] 1.0000000 0.5000000 0.5000000
[2,] 0.7071068 0.7071068 0.4000000
[3,] 0.6324555 0.7905694 0.5757359
[4,] 0.7587726 0.6589589 0.6205267
[5,] 0.7877352 0.6347311 0.5447364
[6,] 0.7380626 0.6774493 0.5273589
[7,] 0.7261948 0.6885205 0.5571624
[8,] 0.7464331 0.6698524 0.5642831
[9,] 0.7511878 0.6656125 0.5521401
[10,] 0.7430613 0.6728920 0.5492873
[11,] 0.7411392 0.6746371 0.5541632
[12,] 0.7444214 0.6716626 0.5553165
[13,] 0.7451956 0.6709648 0.5533472
[14,] 0.7438731 0.6721577 0.5528826
[15,] 0.7435608 0.6724400 0.5536762It seems we are close to convergence after 15 rounds. We are interested in output, Y (the terminal value after convergence).
We want to target a given level \(Y^*=2\) (for instance) that matches, say, observations in the data. But to do so, we only need one parameter. We we wanted to capture more features (inflation, unemployment, etc.), we could potentially use more parameters (if they have a significant impact on the desired outcome).
Below, we allow the function to take \(A\) in vector form (it’s linear basically).
model <- function(A, a, b1, K, nb_sim){
Y <- matrix(1, nrow = nb_sim, ncol = length(A))
w <- matrix(0.5, nrow = nb_sim, ncol = length(A))
N <- matrix(0.5, nrow = nb_sim, ncol = length(A))
for (t in 2:nb_sim){
Y[t,] <- A * (K^a) * N[t-1,]^(1-a) # (1) Cobb-Douglass production function
w[t,] <- A * (1-a) * (K^a) * N[t-1,]^(-a) # (2) Labor demand
N[t,] <- 1 - b1 / w[t-1,] # (3) Labor supply
}
Y[nb_sim,]
}
model(A = 1:5, a = 0.5, b1 = 0.3, K = 1, nb_sim = 20)[1] 0.744041 1.722375 2.714963 3.711234 4.708992ggplot() + xlim(1, 5) + geom_function(fun = model,
args = list(a = 0.5, b1 = 0.3, K = 1, nb_sim = 15)) +
theme_classic() +
geom_hline(yintercept = 2, linetype = 2, color = "gray")
If we have a target production of \(Y^*\), then the value of \(A\) that satisfies this condition is where the two lines meet. Note that in some case, it may be possible that \(\boldsymbol{\theta}\) does not allow to perfectly match the data (desired targets).
In this case, one must resort to error minimization:
\[\underset{\boldsymbol{\theta}}{\min} \ \sum_\textbf{S} || f(\textbf{S, }\boldsymbol{\theta}) - f_{\text{observed}}(\textbf{S}) ||. \]
What happens if we have 2 degrees of freedom? Say \(A\) and \(a\)… We need to plot the outcome in 3D.
model <- function(A, a, b1, K, nb_sim){
Y <- matrix(1, nrow = nb_sim, ncol = length(A))
w <- matrix(0.5, nrow = nb_sim, ncol = length(A))
N <- matrix(0.5, nrow = nb_sim, ncol = length(A))
for (t in 2:nb_sim){
Y[t,] <- A * (K^a) * N[t-1,]^(1-a) # (1) Cobb-Douglass production function
w[t,] <- A * (1-a) * (K^a) * N[t-1,]^(-a) # (2) Labor demand
N[t,] <- 1 - b1 / w[t-1,] # (3) Labor supply
}
Y[nb_sim,]
}
A <- seq(0.5, 2, by = 0.1)
a <- seq(0.1, 0.9, by = 0.02)
pars <- expand_grid(A, a)
AA <- pars$A
aa <- pars$a
prod <- model(AA, aa, 0.3, 1, 15) |> matrix(ncol = length(A))
plot_ly(type = 'surface',
x = ~A,
y = ~a,
z = prod)Here, we have only focused on production, but we could plot \(w\) as well… And in 2D it’s also possible!
data.frame(prod = model(AA, aa, 0.3, 1, 15)) |>
bind_cols(pars) |>
ggplot(aes(x = a, y = prod, color = as.factor(A))) + geom_line() +
theme_classic() + labs(color = "A") +
theme(axis.title.y = element_blank())
It shows that some parameter values do not lead to convergence…