24  Calibrating a Model with the Method of Simulated Moments

Authors

Karsten Kohler

Nitin Nair

Overview

So far, all models on this website have been parameterised by hand without dedicating much attention to empirical accuracy. For example, in the IS-LM model in Chapter 4, we may have set the marginal propensity to consume in the Keynesian consumption function (Equation 4.2) to a value of \(c_1=0.6\) without checking that this value actually generates empirically plausible results. This raises the question of how we can set a model’s parameter values so that its simulated data match empirical data.

This chapter presents a technique to calibrate a model’s parameters against empirical data called the method of simulated moments (MSM) (Duffie and Singleton 1993; McFadden 1989). The MSM chooses parameters to minimise the distance between a set of data moments (e.g. means, variances, autocovariances, cross-covariances, …) and the same moments computed from the model’s simulated data. To stick with the Keynesian consumption function example, we might use empirical data on the correlation between aggregate output and consumption for an economy of interest to calibrate the marginal propensity to consume, so that the correlation between aggregate output and consumption in the simulated data of the IS-LM model is as close as possible to the empirically observed one.

The MSM is the workhorse estimation method in the agent-based and behavioural macroeconomics literature, which we follow here.1 We illustrate its application with the New Keynesian 3-equation model covered in Chapter 19.2

The method of simulated moments

Let \(m^{emp}\) be a \((n_m \times 1)\) column vector of \(n_m\) moments computed from empirical data, and let \(m(\theta)\) be the corresponding vector of simulated moments computed from data simulated by the model at parameter vector \(\theta\). Let \(J(\theta)\) be a quadratic loss function that computes the (weighted) sum of squared deviations of the simulated from the empirical moments:
\[ J(\theta) = \big[m(\theta)-m^{emp}\big]' \, W \, \big[m(\theta)-m^{emp}\big], \tag{24.1}\]

where \(W\) is a \((n_m \times n_m)\) weighting matrix that assigns a weight to each squared deviation in the sum.3

The MSM chooses the parameter vector \(\theta\) so as to minimise the value of the loss function: \[ \hat\theta = \arg\min_\theta \; J(\theta). \] The MSM is typically applied to models that do not have a closed-form solution for \(m(\theta)\), the function that describes how the parameters determine the simulated moments. As a result, \(J(\theta)\) does not have a closed-form solution either and needs to be minimised numerically. In practice, this is done by simulating the model many times for different values of the parameter vector \(\theta\), evaluating the loss function for each of those parameterisations, and eventually choosing the parameterisation that gives the lowest value of the loss function. To smooth out noise, the simulated moments are typically computed as averages over multiple Monte Carlo runs of the model for a given parameterisation.

Typically, the theoretically admissible parameter space (i.e. the total number of conceivable parameterisations) is vast – too large to be explored exhaustively. In practice, only a subset of the parameter space will be used to evaluate the loss function. We will discuss below some strategies for exploring the parameter space.

A further practical issue is the choice of which empirical moments to target. Table 1 summarises some of the most commonly used ones, the feature of the data they capture, and when it is useful to include them in the vector of target moments \(m^{emp}\).

Table 1: Common target moments in the MSM

Moment What it captures When is it useful
Mean Average levels To match the observed level of an individual variable
Variance Volatility To match the observed fluctuations of an individual variable
Autocovariance Persistence To match how long an observed change in an individual variable lasts
Cross-covariance Comovement To match the observed contemporaneous relationship between two variables
Lagged cross-covariance Dynamic relationships To match observed lead-lag relationships between two variables

We will illustrate the choice of moments in the applied example below.

A final practical issue concerns the choice of the weighting matrix \(W\). Three common choices, in increasing order of sophistication, are:

  • Identity weighting \((W=I)\): treats every moment as equally important, but ignores the fact that moments are typically measured with very different precision, may be measured in different units or are naturally at different scales.
  • Diagonal inverse-variance weighting: use a diagonal \(W\) with entries \(1/\widehat{\mathrm{Var}}(m^{emp}_i)\), so that imprecisely estimated moments or moments at naturally larger scales than others are down-weighted.
  • Fully efficient weighting \((W=\hat\Omega^{-1})\): uses the inverse of the entire variance-covariance matrix \(\Omega\) of the empirical moments, which is the asymptotically optimal weighting. It takes into account not only the variances of the moments, but also their covariances.

We will discuss these different options for the weighting matrix in the applied example below.

Applied example: calibrating the New Keynesian 3-equation model

Recall the New Keynesian 3-equation model from Carlin and Soskice (2014), presented in Chapter 19:

\[y_t = A - a_1 r_{t-1} \qquad \text{(IS curve)}\] \[\pi_t = \pi_{t-1} + a_2(y_t - y_e) \qquad \text{(Phillips curve)}\] \[r_t = r_s + a_3(\pi_t - \pi^T) \qquad \text{(Monetary policy rule)}\]

\[r_s = (A-y_e)/a_1 \qquad \text{(stabilising interest rate)} \]

\[a_3 = \frac{1}{a_1\left(\frac{1}{a_2 b} + a_2\right)} \qquad \text{(monetary policy parameter)},\] where \(y\), \(A\), \(r\), \(\pi\), \(y_e\), \(r_s\), and \(\pi^T\) are real output, autonomous demand, the real interest rate, inflation, equilibrium output, the stabilising real interest rate, and the inflation target, respectively. The parameters \(a_1, a_2\) and \(b\) are the sensitivity of aggregate demand with respect to the interest rate, the sensitivity of inflation to the output gap, and the sensitivity of the central bank to deviations of inflation from target.

This version of the model is deterministic: it does not contain any random shocks. To calibrate it against noisy real-world macro data, we need to (i) decide which of its parameters are free to calibrate, and (ii) add a stochastic structure so the model actually produces simulated variances and covariances we can compare against data.

Choosing the parameters to be calibrated

A first important question when calibrating a model is which parameters are to be calibrated and which are fixed externally. The 3-equation model has six exogenous parameters/variables: \((A, y_e, \pi^T, a_1,a_2,b)\). We can fix the inflation target externally at a value of 2% followed explicitly or implicitly by many central banks, i.e. \(\pi^T=2\). Setting autonomous demand \(A\) and potential output \(y_e\) to realistic values is less straightforward as these are not directly observable. Suppose we are not interested in reaching empirically realistic average values for output and the real interest rate in levels. Let us thus fix \(A = y_e = 0\), so that \(r_s=0\), and the steady state becomes \((y^*,\pi^*,r^*) = (0, \pi^T, 0)\). This means \(y_t\) and \(r_t\) can be interpreted as gaps with clean empirical counterparts (an output gap and a de-meaned real interest rate). This normalisation implies that the first moments (= means) of \(y,\pi,r\) carry no information about the remaining three free parameters \(a_1,a_2,b\) and should not be used as target moments. Instead, we will only use variances, autocovariances, and cross-covariances to identify them.

Next, we need to introduce stochasticity to the model. To this end, three shocks are added, one per equation – a demand shock in the IS curve \((\varepsilon_{d,t})\), a cost-push shock in the Phillips curve \((\varepsilon_{s,t})\), and a monetary-policy shock in the policy rule \((\varepsilon_{m,t})\). Along with the aforementioned parameter restrictions, this gives the following stochastic version of the 3-equation model:

\[y_t = -a_1 r_{t-1} + \varepsilon_{d,t}\] \[\varepsilon_{d,t} = \rho_d \varepsilon_{d,t-1} + \mu_{d,t}\] \[\pi_t = \pi_{t-1} + a_2 y_t + \varepsilon_{s,t}\] \[r_t = a_3(\pi_t - \pi^T) + \varepsilon_{m,t},\]

where we let the demand shock \(\varepsilon_{d,t}\) follow an AR(1) process so as to introduce some persistence, which is a common property of macroeconomic time series. The other shocks are assumed to be i.i.d. for simplicity and to save on parameters. All random terms are normally distributed with mean zero: \(\mu_{d,t} \sim N(0,\sigma_d^2)\), \(\varepsilon_{s,t} \sim N(0,\sigma_s^2)\), and \(\varepsilon_{m,t} \sim N(0,\sigma_m^2)\). The shock processes thus contain the parameters \(\rho_d, \sigma_d, \sigma_s,\) and \(\sigma_m\), which are the autocorrelation coefficient of the demand shock, and the standard deviations of the demand, cost-push, and monetary policy shocks, respectively.

This gives \(n_p=7\) free parameters to be calibrated: \(\theta = (a_1, a_2, b, \rho_d, \sigma_d, \sigma_s, \sigma_m)\), with stability requiring \(a_1,a_2,b>0\) and \(|\rho_d|<1\).

Data and empirical moments

The next important step is to choose an appropriate set of empirical moments. The number of moments should generally be at least as large as the number of parameters \((n_m \geq n_p)\). The key variables in the 3-equation model are the output gap, the inflation rate, and the policy rate. We will use their variances, first-order autocovariances, and contemporaneous cross-covariances, giving 9 target moments to calibrate 7 parameters. We will use data for the US economy.

To construct an empirical measure of the output gap, we use the Hamilton regression filter (Hamilton 2018).4 It regresses \(y_{t+h}\) on a constant and \(p\) lags \(y_t,\dots,y_{t-p+1}\), yielding regression residuals that represent the cyclical component. We use the standard settings for quarterly data of \(h=8\) and \(p=4\).

Inflation is annualised quarter-on-quarter core PCE inflation, and the real policy rate is the quarterly-averaged effective federal funds rate minus inflation, de-meaned so its sample mean is exactly zero (matching the normalised steady state of \(r^*=0\)). We restrict the sample to 1996Q1-2007Q4 to focus on the Fed’s de facto 2% inflation-targeting era and to exclude the post-2008 zero-lower bound period, where the policy rate was not really reacting to inflation in the sense the model assumes.

The first block of code loads the relevant US time series data from FRED and constructs the desired series (output gap, inflation rate, and de-meaned real policy rate). Note that FRED data are continuously extended and occasionally revised. Since the Hamilton filter is estimated on the full available sample, the resulting output gap and its moments may change slightly whenever new data come in. To keep the results in this chapter reproducible, the code reads a snapshot of the raw series that was downloaded from FRED on 13 September 2026.5 Setting use_snapshot = FALSE downloads the latest data vintage from FRED instead.

# Clear the environment
rm(list = ls())

# Load relevant libraries
library(quantmod)
library(xts)
library(zoo)

# Suppress an unnecessary warning
options(xts.message.period.apply.mean = FALSE)

# Set desired sample period
sample_start = as.yearqtr("1996 Q1")
sample_end   = as.yearqtr("2007 Q4")

# Choose whether to use the stored data snapshot (which reproduces the results in this
# chapter) or to download the latest data vintage from FRED 
use_snapshot = TRUE

if (use_snapshot) {
  
  # Read snapshot of the raw series (downloaded from FRED on 13 September 2026), either
  # from a local copy or, if there is none, from the website's repository
  snapshot = "data/msm_us_data_raw.csv"
  if (!file.exists(snapshot)) {
    snapshot = paste0("https://raw.githubusercontent.com/DIY-Macro-Sim/",
                      "DIY-Macro-sim-website/main/", snapshot)
  }
  data_raw = read.csv(snapshot)
  data_raw$qtr = as.yearqtr(data_raw$qtr)
  
} else {
  
  # Pull relevant series from FRED (real GDP, core PCE price index, nominal federal funds rate)
  getSymbols(c("GDPC1", "PCEPILFE", "FEDFUNDS"), src = "FRED")

  # Align the three series to quarterly frequency and drop missing values
  gdp      = na.omit(GDPC1)                       
  pcepilfe = apply.quarterly(na.omit(PCEPILFE), mean)  # monthly -> quarterly average
  fedfunds = apply.quarterly(na.omit(FEDFUNDS), mean)  # monthly -> quarterly average

  # Put together data frame with the raw series
  data_raw = data.frame(qtr = as.yearqtr(index(gdp)), gdp = as.numeric(gdp))
  data_raw = merge(data_raw, data.frame(qtr = as.yearqtr(index(pcepilfe)),
                                         pcepilfe = as.numeric(pcepilfe)), by = "qtr", all = TRUE)
  data_raw = merge(data_raw, data.frame(qtr = as.yearqtr(index(fedfunds)),
                                         fedfunds = as.numeric(fedfunds)), by = "qtr", all = TRUE)
  data_raw = data_raw[order(data_raw$qtr), ]
}

# Write function to apply the Hamilton filter and return the cyclical component
# NB: need to check which possible regression observations have no missing or infinite values
hamilton_filter = function(y, h = 8, p = 4) {
  n = length(y)
  t_all = p:(n - h)
  valid = sapply(t_all, function(tt) all(is.finite(c(y[tt + h], y[(tt - p + 1):tt]))))
  t = t_all[valid]
  X = cbind(1, sapply(0:(p - 1), function(j) y[t - j]))
  model = lm(y[t + h] ~ X - 1)
  cycle = rep(NA_real_, n)
  cycle[t + h] = residuals(model)
  cycle
}

# Put together data frame for the transformed series
data_transformed = data.frame(qtr = data_raw$qtr)
data_transformed$y  = 100 * hamilton_filter(log(data_raw$gdp))    # output gap
data_transformed$pi = 400 * c(NA, diff(log(data_raw$pcepilfe)))   # annualised inflation
data_transformed$r  = data_raw$fedfunds - data_transformed$pi     # real policy rate

# Restrict to the desired sample period and remove missing values
data_transformed = data_transformed[data_transformed$qtr >= sample_start &
                                     data_transformed$qtr <= sample_end, ]
data_transformed = na.omit(data_transformed)

# Demean the real policy rate
data_transformed$r = data_transformed$r - mean(data_transformed$r)

Note that the Python code will yield somewhat different numerical results from the R code whose output is presented and discussed in the text. The two languages use different random number generators and different implementations of the minimisation algorithm, and the calibration turns out to be sensitive to such differences.

# Load relevant libraries
import os
import numpy as np
import pandas as pd

# Set desired sample period
sample_start = pd.Period("1996Q1", freq="Q")
sample_end   = pd.Period("2007Q4", freq="Q")

# Choose whether to use the stored data snapshot or to download the latest data
# vintage from FRED
use_snapshot = True

if use_snapshot:

    # Read snapshot of the raw series (downloaded from FRED on 13 September 2026), either
    # from a local copy or, if there is none, from the website's repository
    snapshot = "data/msm_us_data_raw.csv"
    if not os.path.exists(snapshot):
        snapshot = ("https://raw.githubusercontent.com/DIY-Macro-Sim/"
                    "DIY-Macro-sim-website/main/" + snapshot)
    data_raw = pd.read_csv(snapshot)
    data_raw["qtr"] = pd.PeriodIndex(data_raw["qtr"].str.replace(" ", ""), freq="Q")

else:

    # Define function to pull a series from FRED and drop missing values
    def get_fred(series):
        url = "https://fred.stlouisfed.org/graph/fredgraph.csv?id=" + series
        df = pd.read_csv(url, index_col=0, parse_dates=True, na_values=".")
        return df[series].dropna()

    # Pull relevant series from FRED (real GDP, core PCE price index, nominal federal
    # funds rate) and align them to quarterly frequency
    gdp      = get_fred("GDPC1")
    pcepilfe = get_fred("PCEPILFE")
    fedfunds = get_fred("FEDFUNDS")
    gdp.index = gdp.index.to_period("Q")
    pcepilfe  = pcepilfe.groupby(pcepilfe.index.to_period("Q")).mean()  # monthly -> quarterly average
    fedfunds  = fedfunds.groupby(fedfunds.index.to_period("Q")).mean()  # monthly -> quarterly average

    # Put together data frame with the raw series
    data_raw = pd.concat([gdp, pcepilfe, fedfunds], axis=1,
                         keys=["gdp", "pcepilfe", "fedfunds"]).sort_index()
    data_raw = data_raw.rename_axis("qtr").reset_index()

# Write function to apply the Hamilton filter and return the cyclical component
# NB: need to check which possible regression observations have no missing or infinite values
def hamilton_filter(y, h=8, p=4):
    y = np.asarray(y, dtype=float)
    n = len(y)
    t_all = np.arange(p - 1, n - h)
    valid = [np.all(np.isfinite(np.r_[y[tt + h], y[tt - p + 1:tt + 1]])) for tt in t_all]
    t = t_all[valid]
    X = np.column_stack([np.ones(len(t))] + [y[t - j] for j in range(p)])
    beta = np.linalg.lstsq(X, y[t + h], rcond=None)[0]
    cycle = np.full(n, np.nan)
    cycle[t + h] = y[t + h] - X @ beta
    return cycle

# Put together data frame for the transformed series
data_transformed = pd.DataFrame({"qtr": data_raw["qtr"]})
data_transformed["y"]  = 100 * hamilton_filter(np.log(data_raw["gdp"]))  # output gap
data_transformed["pi"] = 400 * np.log(data_raw["pcepilfe"]).diff()       # annualised inflation
data_transformed["r"]  = data_raw["fedfunds"] - data_transformed["pi"]   # real policy rate

# Restrict to the desired sample period and remove missing values
data_transformed = data_transformed[(data_transformed["qtr"] >= sample_start) &
                                    (data_transformed["qtr"] <= sample_end)]
data_transformed = data_transformed.dropna()

# Demean the real policy rate
data_transformed["r"] = data_transformed["r"] - data_transformed["r"].mean()

Next, we write a function to compute the nine desired empirical moments from the \(y\), \(\pi\), and \(r\) series:

# Write function that computes the nine desired moments from y, pi, and r series
compute_moments = function(mat) {
  y = mat[, "y"]; pi = mat[, "pi"]; r = mat[, "r"]
  n = nrow(mat)
  c(
    var_y    = var(y),
    var_pi   = var(pi),
    var_r    = var(r),
    ar1_y    = cov(y[2:n],  y[1:(n - 1)]),
    ar1_pi   = cov(pi[2:n], pi[1:(n - 1)]),
    ar1_r    = cov(r[2:n],  r[1:(n - 1)]),
    cov_y_pi = cov(y, pi),
    cov_y_r  = cov(y, r),
    cov_pi_r = cov(pi, r)
  )
}

# Convert the empirical series to a matrix and apply function to compute empirical moments
ts_data = as.matrix(data_transformed[, c("y", "pi", "r")])
obs_moments_raw = compute_moments(ts_data)

# Print rounded moments
round(obs_moments_raw, 2)
   var_y   var_pi    var_r    ar1_y   ar1_pi    ar1_r cov_y_pi  cov_y_r 
    3.94     0.28     3.53     3.71     0.08     3.30    -0.12     2.29 
cov_pi_r 
   -0.31 
# Write function that computes the nine desired moments from a matrix whose columns
# contain the y, pi, and r series
moment_names = ["var_y", "var_pi", "var_r", "ar1_y", "ar1_pi", "ar1_r",
                "cov_y_pi", "cov_y_r", "cov_pi_r"]

def cov(x, z):
    return np.cov(x, z)[0, 1]

def compute_moments(mat):
    y = mat[:, 0]; pi = mat[:, 1]; r = mat[:, 2]
    return np.array([
        np.var(y, ddof=1),
        np.var(pi, ddof=1),
        np.var(r, ddof=1),
        cov(y[1:], y[:-1]),
        cov(pi[1:], pi[:-1]),
        cov(r[1:], r[:-1]),
        cov(y, pi),
        cov(y, r),
        cov(pi, r)
    ])

# Convert the empirical series to a matrix and apply function to compute empirical moments
ts_data = data_transformed[["y", "pi", "r"]].to_numpy()
obs_moments_raw = compute_moments(ts_data)

# Print rounded moments
print(pd.Series(obs_moments_raw, index=moment_names).round(2))

Simulated moments

Next, we write a function that simulates the stochastic version of the model. The function takes as arguments the parameter vector \((\theta)\), the desired time horizon of the simulation (T_sim), and a chosen “burn-in” period of initial data points that will be discarded (burnin). Three further arguments specify the implementation of a manual shock scenario that will be used later to create illustrative impulse response functions. With the default setting shock_t = NULL, we can ignore these arguments for now.

# Set inflation target externally
pt = 2  

# Define function to simulate the model
simulate_model = function(theta, T_sim = 300, burnin = 100,
                            shock_t = NULL, shock_size = 0, shock_eq = "d") {
  
  # Assign parameter values from parameter vector to individual parameters
  a1 = theta[1] 
  a2 = theta[2] 
  b = theta[3]
  rho_d = theta[4]
  sigma_d = theta[5] 
  sigma_s = theta[6] 
  sigma_m = theta[7]
  a3 = (a1 * (1 / (a2 * b) + a2))^(-1)

  # Create vectors that will contain the simulated data
  y = numeric(T_sim); p = numeric(T_sim); r = numeric(T_sim); eps_d = numeric(T_sim)
  
  # Create shock processes
  mu_d = rnorm(T_sim, 0, sigma_d)
  shock_s = rnorm(T_sim, 0, sigma_s)
  shock_m = rnorm(T_sim, 0, sigma_m)

  # Create (optional) additional shock scenario (will later be used to generate IRF)
  if (!is.null(shock_t)) {
    if (shock_eq == "d") mu_d[shock_t] = mu_d[shock_t] + shock_size
    if (shock_eq == "s") shock_s[shock_t] = shock_s[shock_t] + shock_size
    if (shock_eq == "m") shock_m[shock_t] = shock_m[shock_t] + shock_size
  }

  # Initialise at steady state
  y[1] = 0; p[1] = pt; r[1] = 0

  # Simulate
  for (t in 2:T_sim) {
    eps_d[t] = rho_d * eps_d[t-1] + mu_d[t]       # AR(1) demand shock
    y[t] = -a1 * r[t - 1] + eps_d[t]            # IS curve, A=0
    p[t] = p[t - 1] + a2 * y[t] + shock_s[t]    # Phillips curve, ye=0
    r[t] = a3 * (p[t] - pt) + shock_m[t]        # Policy rule, rs=0
  }

  # Return simulated data for the chosen simulation period
  res = cbind(y = y, p = p, r = r)
  return(res[(burnin + 1):T_sim, ])
}
# Set inflation target externally
pt = 2

# Define function to simulate the model
# NB: 'rng' is the random number generator from which the shocks are drawn
def simulate_model(theta, T_sim=300, burnin=100, rng=None,
                   shock_t=None, shock_size=0, shock_eq="d"):

    # Assign parameter values from parameter vector to individual parameters
    a1, a2, b, rho_d, sigma_d, sigma_s, sigma_m = theta
    a3 = (a1 * (1 / (a2 * b) + a2))**(-1)

    # Create vectors that will contain the simulated data
    y = np.zeros(T_sim); p = np.zeros(T_sim); r = np.zeros(T_sim); eps_d = np.zeros(T_sim)

    # Create shock processes
    if rng is None:
        rng = np.random.default_rng()
    mu_d    = rng.normal(0, sigma_d, T_sim)
    shock_s = rng.normal(0, sigma_s, T_sim)
    shock_m = rng.normal(0, sigma_m, T_sim)

    # Create (optional) additional shock scenario (will later be used to generate IRF)
    if shock_t is not None:
        if shock_eq == "d": mu_d[shock_t] += shock_size
        if shock_eq == "s": shock_s[shock_t] += shock_size
        if shock_eq == "m": shock_m[shock_t] += shock_size

    # Initialise at steady state
    y[0] = 0; p[0] = pt; r[0] = 0

    # Simulate
    for t in range(1, T_sim):
        eps_d[t] = rho_d * eps_d[t - 1] + mu_d[t]   # AR(1) demand shock
        y[t] = -a1 * r[t - 1] + eps_d[t]            # IS curve, A=0
        p[t] = p[t - 1] + a2 * y[t] + shock_s[t]    # Phillips curve, ye=0
        r[t] = a3 * (p[t] - pt) + shock_m[t]        # Policy rule, rs=0

    # Return simulated data (columns: y, p, r) for the chosen simulation period
    res = np.column_stack([y, p, r])
    return res[burnin:, :]

Then we define another function that computes the simulated moments for a given parameterisation as averages over \(MC\) Monte Carlo runs (which we set to \(MC=50\) here). By averaging the moments over multiple runs, we reduce sensitivity of the results to outlier values. Note that we fix a random seed inside the sim_moments() function, so it resets to the same reproducible sequence of random numbers every time the function is called - which will happen repeatedly below when we simulate the model for different parameterisations. This means that while every individual Monte Carlo run performed via the simulate_model function is based on different random values for the shocks, the sequence of Monte Carlo runs performed for every parameterisation uses the same underlying sequence of shocks; only the parameterisation changes what happens to them. This approach keeps the loss function comparable and smooth across evaluations.

# Write function that computes simulated moments as Monte Carlo averages
sim_moments = function(theta, MC = 50) {
  
  # fix random seed
  set.seed(123)                 
  
  # create matrix in which simulated moments will be stored
  all_moms = matrix(NA, MC, 9)  
  
  for (s in 1:MC) {             # Monte Carlo loop
    sim_data = simulate_model(theta)
    y_s = sim_data[, "y"]; p_s = sim_data[, "p"]; r_s = sim_data[, "r"]
    n = nrow(sim_data)
    
    # Compute simulated moments  
    all_moms[s, ] = c(
      var(y_s), var(p_s), var(r_s),
      cov(y_s[2:n], y_s[1:(n - 1)]),
      cov(p_s[2:n], p_s[1:(n - 1)]),
      cov(r_s[2:n], r_s[1:(n - 1)]),
      cov(y_s, p_s), cov(y_s, r_s), cov(p_s, r_s)
    )
   }                          # close loop
  
  # Return moments
  moms = colMeans(all_moms, na.rm = TRUE)
  if (any(!is.finite(moms))) return(rep(1e6, 9))
  return(moms)
}
# Write function that computes simulated moments as Monte Carlo averages
def sim_moments(theta, MC=50):

    # fix random seed
    rng = np.random.default_rng(123)

    # create matrix in which simulated moments will be stored
    all_moms = np.full((MC, 9), np.nan)

    for s in range(MC):           # Monte Carlo loop
        sim_data = simulate_model(theta, rng=rng)

        # Compute simulated moments (re-using the function defined above)
        all_moms[s, :] = compute_moments(sim_data)

    # Return moments
    moms = np.nanmean(all_moms, axis=0)
    if np.any(~np.isfinite(moms)):
        return np.full(9, 1e6)
    return moms

Loss function

Then we define a function that computes the value of the loss function (Equation 24.1) for a given set of parameters \(\theta\) and a weighting matrix \(W\) that need to be supplied as arguments to the function:

# Define function that computes the value of the loss function
objective = function(theta, Wmat) {
  sim_moms = sim_moments(theta)
  if (any(!is.finite(sim_moms))) return(1e10)
  diff = sim_moms - obs_moments_raw
  as.numeric(t(diff) %*% Wmat %*% diff)
}
# Define function that computes the value of the loss function
def objective(theta, Wmat):
    sim_moms = sim_moments(theta)
    if np.any(~np.isfinite(sim_moms)):
        return 1e10
    diff = sim_moms - obs_moments_raw
    return float(diff @ Wmat @ diff)

Exploring the parameter space

A key challenge when using a numerical solution method like the MSM is to decide on a procedure for how to explore the, typically vast, parameter space. We will tackle this problem in two ways. First, we utilise R’s general-purpose numerical minimiser, optim(). This is not specifically designed for the MSM – it just repeatedly evaluates whatever objective function you supply at different trial parameter vectors, and searches for the one that returns the smallest number. Here we use method="L-BFGS-B", a quasi-Newton algorithm: at each trial point it numerically approximates the gradient — the local slope of the objective function, i.e. by how much the loss would change if each parameter were nudged up or down slightly, and in which direction that loss decreases fastest. It then uses accumulated gradient information to approximate the objective function’s local curvature to decide which parameterisation to evaluate next. The optim() function takes as arguments two constraints (lower/upper) that specify the admissible min and max values for each parameter. Specifying these bounds also allows us to enforce \(a_1,a_2,b>0\) and \(|\rho_d|<1\), which we require to ensure stability of the stochastic version of the 3-equation model. Finally, the control argument determines when the algorithm stops: after at most maxit iterations, or once the relative improvement in the loss between two iterations falls below factr times machine precision (here about \(2 \times 10^{-4}\)), or once the largest component of the gradient falls below pgtol. We choose relatively loose tolerances to keep the computation time manageable.

Two things are worth flagging about optim(). First, it returns a convergence code: 0 means it believes it stopped at a genuine local minimum, while nonzero codes flag a problem. Second, and more importantly: optim() is a local search. It has no way of knowing whether the minimum it finds is the best one available, or just the closest one to wherever it started. Consequently, the result might be sensitive to the initial parameterisation that needs to be supplied to the function.

This brings us to the second way in which we will tackle the parameter choice problem. In a robustness test, which we perform below, we will not only consider one initial parameterisation, but many that are sampled from across the admissible parameter space. This allows us to assess the sensitivity of the optimal parameterisation to the initial values, and to potentially find better parameterisations.

Taking a first shot at the calibration

For our first shot at the calibration, we will use a simplified approach. First, we will start the minimisation algorithm from only a single initial parameterisation. Second, we will use identity weighting in the loss function, i.e. we will set \(W=I\). This will produce a tentative ‘optimal’ parameterisation. We will then refine our approach to see if it yields a different, potentially better, parameterisation.

# Define an initial parameterisation
theta_start = c(a1 = 0.3, a2 = 0.7, b = 1, rho_d = 0.5,
                 sigma_d = 0.3, sigma_s = 0.3, sigma_m = 0.3)

# Define lower and upper bounds for the admissible parameter values
bounds_lower = c(0.01, 0.01, 0.01, -0.99, 0.001, 0.001, 0.001)
bounds_upper = c(2, 2, 5, 0.99, 5, 5, 5)

# Set up an identity weighting matrix
W_identity = diag(length(obs_moments_raw))    

# Run minimisation algorithm on loss function
fit_id = optim(theta_start, function(th) objective(th, W_identity),
                method = "L-BFGS-B",
                lower = bounds_lower, upper = bounds_upper,
                control = list(maxit = 1000, factr = 1e12, pgtol = 1e-3))

# Print fit
cat("Convergence:", fit_id$convergence, " Objective:", round(fit_id$value, 4), "\n")
Convergence: 0  Objective: 7.398 
# Print calibrated parameters (including the implied a3 parameter)
theta_id = fit_id$par
a3_id = with(as.list(theta_id), (a1 * (1 / (a2 * b) + a2))^(-1))
round(c(theta_id, a3 = a3_id), 4)
     a1      a2       b   rho_d sigma_d sigma_s sigma_m      a3 
 0.0100  0.0203  2.0460  0.5367  1.8726  0.0115  0.5202  4.1568 
from scipy.optimize import minimize

# Define an initial parameterisation
theta_names = ["a1", "a2", "b", "rho_d", "sigma_d", "sigma_s", "sigma_m"]
theta_start = np.array([0.3, 0.7, 1, 0.5, 0.3, 0.3, 0.3])

# Define lower and upper bounds for the admissible parameter values
bounds_lower = np.array([0.01, 0.01, 0.01, -0.99, 0.001, 0.001, 0.001])
bounds_upper = np.array([2, 2, 5, 0.99, 5, 5, 5])
bounds = list(zip(bounds_lower, bounds_upper))

# Set up an identity weighting matrix
W_identity = np.eye(len(obs_moments_raw))

# Set options for the minimisation algorithm (mirroring those used in the R code)
opts = {"maxiter": 1000, "ftol": 1e12 * np.finfo(float).eps, "gtol": 1e-3, "eps": 1e-3}

# Run minimisation algorithm on loss function
fit_id = minimize(objective, theta_start, args=(W_identity,),
                  method="L-BFGS-B", bounds=bounds, options=opts)

# Print fit
print("Convergence:", fit_id.status, " Objective:", round(fit_id.fun, 4))

# Define function to compute the implied a3 parameter
def compute_a3(theta):
    a1, a2, b = theta[0], theta[1], theta[2]
    return (a1 * (1 / (a2 * b) + a2))**(-1)

# Print calibrated parameters (including the implied a3 parameter)
theta_id = fit_id.x
a3_id = compute_a3(theta_id)
print(pd.Series(np.append(theta_id, a3_id), index=theta_names + ["a3"]).round(4))

Before evaluating the goodness of fit of this tentative parameterisation, let us check if we’d get a different result when using a different weighting matrix.

Taking a second shot with a diagonal weighting matrix

Next, we re-run the calibration using a more sophisticated version of the weighting matrix \(W\). To weight moments by how precisely they are estimated, we need to estimate the variance-covariance matrix of the empirical moments \((\Omega)\). We estimate \(\Omega\) with a block bootstrap of the empirical data, using the tsboot() function from R’s boot package. Rather than resampling individual quarters independently, a block bootstrap resamples (with replacement) contiguous blocks of quarters, which preserves the autocorrelation in the data. We use the specification sim="geom", which draws blocks of random (geometrically distributed) length rather than a single fixed length, with a mean block length of l=8 quarters (= two years). We perform R=5000 bootstrap replications.

Note that in order to call the tsboot() function, we need to supply a function as an argument that computes the statistics that are to be bootstrapped. For this, we re-use the previously created compute_moments function.

library(boot)

# Set random seed for reproducibility
set.seed(123)

# Create a bootstrap sample of the moments
boot_out  = tsboot(ts_data, compute_moments, R = 5000, l = 8, sim = "geom")

# Save bootstrap sample as a matrix
boot_moms = boot_out$t
colnames(boot_moms) = names(obs_moments_raw)

# Compute variance-covariance matrix of the bootstrapped moments sample
Omega_boot = cov(boot_moms)

# Print estimated standard errors of the moments
round(sqrt(diag(Omega_boot)), 2)
   var_y   var_pi    var_r    ar1_y   ar1_pi    ar1_r cov_y_pi  cov_y_r 
    1.25     0.04     0.92     1.15     0.04     0.90     0.20     1.03 
cov_pi_r 
    0.14 
# Write function that performs a block bootstrap with blocks of random (geometrically
# distributed) length, where 'l' is the mean block length and 'R' the number of replications
def block_bootstrap(data, statistic, R=5000, l=8, seed=123):

    # Set random seed for reproducibility
    rng = np.random.default_rng(seed)

    n = data.shape[0]
    boot_stats = np.empty((R, len(statistic(data))))

    for i in range(R):
        idx = np.empty(n, dtype=int)
        filled = 0
        while filled < n:
            start  = rng.integers(n)                 # random starting point of the block
            length = rng.geometric(1 / l)            # random length of the block
            block  = (start + np.arange(length)) % n # wrap around at the end of the sample
            m = min(length, n - filled)
            idx[filled:filled + m] = block[:m]
            filled += m
        boot_stats[i, :] = statistic(data[idx, :])   # compute statistics on resampled data

    return boot_stats

# Create a bootstrap sample of the moments
boot_moms = block_bootstrap(ts_data, compute_moments, R=5000, l=8, seed=123)

# Compute variance-covariance matrix of the bootstrapped moments sample
Omega_boot = np.cov(boot_moms, rowvar=False)

# Print estimated standard errors of the moments
print(pd.Series(np.sqrt(np.diag(Omega_boot)), index=moment_names).round(2))

Having obtained an estimate for the variance-covariance matrix \(\hat{\Omega}\) of the moments, we re-run the calibration using diagonal inverse-variance weighting, where \(W\) is a diagonal matrix with entries \(1/\widehat{\mathrm{Var}}(m^{emp}_i)\), i.e. the inverse of the main diagonal of \(\hat{\Omega}\).6 Note that we use the previously calibrated optimal parameter vector as an initialisation for the minimisation algorithm.

# Construct diagonal weighting matrix
W_opt_diag = diag(1 / diag(Omega_boot))

# Run minimisation algorithm
fit_diag = optim(theta_id, function(th) objective(th, W_opt_diag),
                   method = "L-BFGS-B", lower = bounds_lower, upper = bounds_upper,
                   control = list(maxit = 1000, factr = 1e12, pgtol = 1e-3))
theta_diag = fit_diag$par

# Compute the implied monetary policy parameter a3 
a3_diag = with(as.list(theta_diag), (a1 * (1 / (a2 * b) + a2))^(-1))

# Print both parameterisations
round(rbind(identity = c(theta_id, a3 = a3_id), diagonal = c(theta_diag, a3 = a3_diag)), 4)
             a1     a2      b  rho_d sigma_d sigma_s sigma_m     a3
identity 0.0100 0.0203 2.0460 0.5367  1.8726  0.0115  0.5202 4.1568
diagonal 0.2938 0.0100 1.4347 0.8506  1.0728  0.0323  1.8468 0.0488
# Construct diagonal weighting matrix
W_opt_diag = np.diag(1 / np.diag(Omega_boot))

# Run minimisation algorithm
fit_diag = minimize(objective, theta_id, args=(W_opt_diag,),
                    method="L-BFGS-B", bounds=bounds, options=opts)
theta_diag = fit_diag.x

# Compute the implied monetary policy parameter a3
a3_diag = compute_a3(theta_diag)

# Print both parameterisations
print(pd.DataFrame([np.append(theta_id, a3_id), np.append(theta_diag, a3_diag)],
                   index=["identity", "diagonal"], columns=theta_names + ["a3"]).round(4))

It can be seen that the diagonal weighting matrix does give a different optimal parameter vector. However, this difference might also stem from having used a different initialisation from which the algorithm may have converged to a different local minimum. As a final exercise, before we evaluate the fit of the different candidate parameterisations, we will run the minimisation algorithm for different initialisations.

Taking a third shot with multiple initialisations

Recall that optim() only explores the basin closest to the initial parameterisation that is being supplied. To check whether that basin is a genuine (near-)global optimum, we draw a space-filling set of starting points with Latin Hypercube Sampling (LHS): stratify each parameter’s range into equal-probability bins, take exactly one random draw per bin, then randomly pair bins across parameters. This covers the 7-dimensional parameter space more evenly than independent uniform draws would, which could easily leave whole regions unsampled by chance. We draw a sample of size 30 with randomLHS() from R’s lhs package, which returns a sample on the unit hypercube \([0,1]^7\). We then use qunif() to rescale each column to the parameter’s chosen bounds. We also append theta_id and theta_diag to the set of starting points (so that the search is guaranteed to yield a parameterisation that is at least as good as the ones we already found). The optim() function is then re-run for every starting sample, keeping the best result. We continue to use the diagonal weighting matrix in the loss function.

library(lhs)

# Set random seed for reproducibility
set.seed(6)

# Set number of parameter samples and save number of parameters
N_starts = 30
n_par    = length(theta_id)

# Draw a Latin Hypercube sample on the unit hypercube
lhs_unit   = randomLHS(N_starts, n_par)

# Rescale each column to the corresponding parameter's chosen bounds defined previously
lhs_starts = lhs_unit
for (j in 1:n_par) {
  lhs_starts[, j] = qunif(lhs_unit[, j], bounds_lower[j], bounds_upper[j])
}
colnames(lhs_starts) = names(theta_id)

# Append the two parameterisations already found (theta_id, theta_diag) as two more
# starting points
all_starts = rbind(lhs_starts, theta_id = theta_id, theta_diag = theta_diag)
N_all = nrow(all_starts)

# Create a matrix that will contain the optimal parameterisation for each start
multistart_par = matrix(NA, N_all, length(theta_id))
colnames(multistart_par) = names(theta_id)

# Create a matrix that will contain the associated value of the loss function
multistart_obj = rep(NA, N_all)

# Calibrate for each start
for (i in 1:N_all) {
  fit_i = tryCatch(
    optim(all_starts[i, ], function(th) objective(th, W_opt_diag),
          method = "L-BFGS-B", lower = bounds_lower, upper = bounds_upper,
          control = list(maxit = 1000, factr = 1e12, pgtol = 1e-3)),
    error = function(e) NULL
  )
  if (!is.null(fit_i)) {
    multistart_par[i, ] = fit_i$par
    multistart_obj[i]   = fit_i$value
  }
}

# Find and print parameterisation that minimises loss function
best_i = which.min(multistart_obj)
round(rbind(single_start = c(theta_diag, obj = fit_diag$value),
            lhs_best     = c(multistart_par[best_i, ], obj = multistart_obj[best_i])), 4)
                 a1   a2      b  rho_d sigma_d sigma_s sigma_m     obj
single_start 0.2938 0.01 1.4347 0.8506  1.0728  0.0323  1.8468 35.9786
lhs_best     0.2937 0.01 1.4346 0.8495  1.0721  0.0323  1.8469 35.9766
from scipy.stats import qmc

# Set number of parameter samples and save number of parameters
N_starts = 30
n_par    = len(theta_id)

# Draw a Latin Hypercube sample on the unit hypercube (with a random seed for reproducibility)
lhs_unit = qmc.LatinHypercube(d=n_par, seed=6).random(N_starts)

# Rescale each column to the corresponding parameter's chosen bounds defined previously
lhs_starts = qmc.scale(lhs_unit, bounds_lower, bounds_upper)

# Append the two parameterisations already found (theta_id, theta_diag) as two more
# starting points
all_starts = np.vstack([lhs_starts, theta_id, theta_diag])
N_all = all_starts.shape[0]

# Create a matrix that will contain the optimal parameterisation for each start
multistart_par = np.full((N_all, n_par), np.nan)

# Create a vector that will contain the associated value of the loss function
multistart_obj = np.full(N_all, np.nan)

# Calibrate for each start
for i in range(N_all):
    try:
        fit_i = minimize(objective, all_starts[i, :], args=(W_opt_diag,),
                         method="L-BFGS-B", bounds=bounds, options=opts)
        multistart_par[i, :] = fit_i.x
        multistart_obj[i]    = fit_i.fun
    except Exception:
        pass

# Find and print parameterisation that minimises loss function
best_i = np.nanargmin(multistart_obj)
print(pd.DataFrame([np.append(theta_diag, fit_diag.fun),
                    np.append(multistart_par[best_i, :], multistart_obj[best_i])],
                   index=["single_start", "lhs_best"], columns=theta_names + ["obj"]).round(4))

It turns out that parameterisation found earlier using a single initialisation had already found the best fit available: the best of the 32 starts (objective \(\approx\) 36) lands at essentially the same point as theta_diag itself, up to negligible numerical noise. Note that this will not always be the case – as discussed above, optim() gives no guarantee of finding a global optimum – so this robustness check is still worth doing.

Comparing the fit of the candidate parameterisations

Since the search across multiple initialisations returned the same parameterisation as the single diagonal-weighting run, we are left with two distinct candidates to compare: identity weighting and diagonal weighting. Let’s compare them by checking how closely each one’s simulated moments match their empirical counterparts. We will assess the empirical fit by (i) checking whether the simulated moment has the same sign as its empirical counterpart, (ii) the percent deviation of the simulated from the empirical moment, and (iii) a moment missing rate (MMR), proposed by Franke (2022), defined as a moment’s deviation from its empirical counterpart as a percentage of the half-width of the moment’s bootstrapped 95% confidence interval.7

The sign match is a useful piece of qualitative information on whether the model gets the direction of a relationship right. While the percent deviation is an intuitive quantitative measure, it may not be directly comparable across moments: it blows up for a moment whose empirical value is close to zero (e.g. cov_y_pi), and it does not distinguish a large miss on a precisely-estimated moment from an equally large miss on a very noisily-estimated one. The MMR addresses these issues. A value of 0 is a perfect match; below 100 means the simulated moment falls inside the empirical 95% confidence interval (obs \(\pm\) 2*SE); above 100 means it falls outside.

# Compute simulated moments under each of the two parameterisations
sim_id   = sim_moments(theta_id)
sim_diag = sim_moments(theta_diag)

# Retrieve bootstrapped standard error of each moment
se_boot = sqrt(diag(Omega_boot))

# Define function to build goodness-of-fit table for each parameterisation
gof_table = function(sim, obs, se) {
  data.frame(
    moment      = names(obs),
    empirical   = round(obs, 3),
    simulated   = round(sim, 3),
    sign_match  = sign(sim) == sign(obs),
    pct_dev     = round(100 * (sim - obs) / obs, 1),
    mmr         = round(100 * abs(sim - obs) / (2 * se), 1),
    row.names = NULL
  )
}

# Call function
table_identity = gof_table(sim_id,   obs_moments_raw, se_boot)
table_diagonal = gof_table(sim_diag, obs_moments_raw, se_boot)
# Compute simulated moments under each of the two parameterisations
sim_id   = sim_moments(theta_id)
sim_diag = sim_moments(theta_diag)

# Retrieve bootstrapped standard error of each moment
se_boot = np.sqrt(np.diag(Omega_boot))

# Define function to build goodness-of-fit table for each parameterisation
def gof_table(sim, obs, se):
    return pd.DataFrame({
        "moment":     moment_names,
        "empirical":  np.round(obs, 3),
        "simulated":  np.round(sim, 3),
        "sign_match": np.sign(sim) == np.sign(obs),
        "pct_dev":    np.round(100 * (sim - obs) / obs, 1),
        "mmr":        np.round(100 * np.abs(sim - obs) / (2 * se), 1)
    })

# Call function and print tables
table_identity = gof_table(sim_id,   obs_moments_raw, se_boot)
table_diagonal = gof_table(sim_diag, obs_moments_raw, se_boot)
print(table_identity)
print(table_diagonal)

Identity weighting

table_identity

Diagonal weighting

table_diagonal

Looking at the fit of the individual moments, we see that variances of \(r\) and \(y\), and the auto-covariance of \(y\) are matched quite well. By contrast, the match for the variance of \(\pi\) is poor, and also for the auto-covariances of \(\pi\) and \(r\), the latter even getting the sign wrong. Recall that we are not using AR(1) process for \(\pi\) and \(r\), which might partly explain this poor fit. The cross-covariances are also matched relatively poorly, with simulated moments being outside of the 95% confidence band and the simulated covariance between \(y\) and \(\pi\) exhibiting the wrong sign.

Averaging the MMR across all 9 moments gives a single overall goodness-of-fit score for each parameterisation,

round(c(identity = mean(table_identity$mmr),
        diagonal = mean(table_diagonal$mmr)), 1)
identity diagonal 
    98.0     79.7 
print(pd.Series({"identity": table_identity["mmr"].mean(),
                 "diagonal": table_diagonal["mmr"].mean()}).round(1))

suggesting that the diagonal weighting matrix overall delivers the better fit.

Overall, this suggests that there is room for improvement when it comes to calibrating the 3-equation model. One might want to consider a different set of moments (e.g. lagged auto-covariances), a different sample period, using different shock processes in the model, or even re-considering the lag structure assumed by the model. We will discuss some more general caveats below.

Analysing a cost-push shock with the calibrated model

A natural use of a calibrated model is to trace out its predicted response to a shock via an impulse response function (IRF). Here we consider a cost-push shock by perturbing the \(\varepsilon_{s,t}\) term in the Phillips curve by one percentage point.8 We then assess how \(\pi_t, r_t, y_t\) evolve relative to their otherwise identical unshocked path.

# Select optimal parameter vector (the best of all starts, which include theta_diag)
theta_eff = multistart_par[best_i, ]

H = 40              # horizon, in quarters after the shock
shock_period = 5    # quarter (within the simulated window) the shock hits

# Use the previously defined simulation function to simulate the model
# once without and once with the shock
set.seed(42); sim_base  = simulate_model(theta_eff, shock_size = 0)
set.seed(42); sim_shock = simulate_model(theta_eff, shock_size = 1,
                                          shock_t = 100 + shock_period, shock_eq = "s")

# Compute impulse response as the deviation of the shocked path from the (otherwise
# identical) baseline path, indexed in quarters since the shock hit
horizon = 0:H
irf_p = sim_shock[shock_period + horizon, "p"] - sim_base[shock_period + horizon, "p"]
irf_r = sim_shock[shock_period + horizon, "r"] - sim_base[shock_period + horizon, "r"]
irf_y = sim_shock[shock_period + horizon, "y"] - sim_base[shock_period + horizon, "y"]

# Plot IRFs
par(mfrow = c(1, 3))
plot(horizon, irf_p, type = "l", lwd = 2, col = "tomato",
     xlab = "Quarters after shock", ylab = "Deviation from baseline (pp)", main = "Inflation")
abline(h = 0, lty = 2, col = "grey50")

plot(horizon, irf_r, type = "l", lwd = 2, col = "tomato",
     xlab = "Quarters after shock", ylab = "Deviation from baseline (pp)", main = "Policy rate")
abline(h = 0, lty = 2, col = "grey50")

plot(horizon, irf_y, type = "l", lwd = 2, col = "tomato",
     xlab = "Quarters after shock", ylab = "Deviation from baseline (pp)", main = "Output gap")
abline(h = 0, lty = 2, col = "grey50")

Impulse responses to a 1-percentage-point inflationary (cost-push) shock
import matplotlib.pyplot as plt

# Select optimal parameter vector (the best of all starts, which include theta_diag)
theta_eff = multistart_par[best_i, :]

H = 40              # horizon, in quarters after the shock
shock_period = 5    # quarter (within the simulated window) the shock hits

# Use the previously defined simulation function to simulate the model
# once without and once with the shock (using the same random seed for both runs)
sim_base  = simulate_model(theta_eff, rng=np.random.default_rng(42), shock_size=0)
sim_shock = simulate_model(theta_eff, rng=np.random.default_rng(42), shock_size=1,
                           shock_t=100 + shock_period, shock_eq="s")

# Compute impulse response as the deviation of the shocked path from the (otherwise
# identical) baseline path, indexed in quarters since the shock hit
horizon = np.arange(H + 1)
irf_y = sim_shock[shock_period + horizon, 0] - sim_base[shock_period + horizon, 0]
irf_p = sim_shock[shock_period + horizon, 1] - sim_base[shock_period + horizon, 1]
irf_r = sim_shock[shock_period + horizon, 2] - sim_base[shock_period + horizon, 2]

# Plot IRFs
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
for ax, irf, title in zip(axes, [irf_p, irf_r, irf_y],
                          ["Inflation", "Policy rate", "Output gap"]):
    ax.plot(horizon, irf, linewidth=2, color="tomato")
    ax.axhline(0, linestyle="--", color="grey")
    ax.set_xlabel("Quarters after shock")
    ax.set_ylabel("Deviation from baseline (pp)")
    ax.set_title(title)
plt.tight_layout()
plt.show()

Inflation jumps by the full size of the shock on impact and then starts to fall, but extremely slowly.9 This is a direct reflection of the calibrated parameter value of \(a_2=0.01\) in the Phillips curve: with the estimated sensitivity of inflation to the output gap being very low, while inflation expectations are strongly adaptive, it takes extremely long for inflation to return to target. The policy rate does move in the right direction – the central bank raises \(r\), which lowers \(y\), which should help pull inflation back down – but both responses are tiny and just as persistent as the inflation shock itself. The soft monetary policy response is due to the relatively low value of \(a_3 = \frac{1}{a_1\left(\frac{1}{a_2 b} + a_2\right)} \approx \frac{1}{0.29\left(\frac{1}{0.01 \times 1.43} + 0.01\right)} \approx 0.05\) in the monetary policy reaction function. While the central bank cares relatively strongly about inflation \((b=1.43)\), it knows that the reaction of inflation to output is weak. As a result, it is not willing to fight inflation more aggressively, as this would require the creation of a very large output gap. Thus, with a weak link between the output gap and inflation, and a relatively weak monetary policy response, even a sustained increase in \(r\) leaves the shock largely unaddressed over any horizon of practical interest.

Caveats

Beyond the specific fit problems already discussed, our exercise illustrates some more general limitations of the MSM worth keeping in mind.

Moments are correlational, not causal. The MSM identifies parameters from unconditional variances and (cross-)covariances. These are descriptive statistical objects that carry no information about a causal direction on their own. By contrast, the parameters that are to be calibrated, e.g. the response of inflation to the output gap, \(a_2>0\), represent theoretically assumed causal effects. An empirically measured moment such as the cross-covariance between \(y\) and \(\pi\) is driven by all sorts of confounding shocks that occurred during the sample period. A positive cost-push shock, for example, may generate a negative correlation between \(y\) and \(\pi\), whereas a positive demand shock may yield a positive one. It thus cannot be taken for granted that an empirically observed correlation is sufficient for pinning down the causal effect of, e.g., \(y\) on \(\pi\), represented by \(a_2\).

The choice of moments is arbitrary. It is a priori unclear which are the best moments to choose to calibrate a given set of parameters. There is no uniquely correct set of moments to target, so any statement about a calibrated parameter is implicitly conditional on the particular moments (and weights) chosen. A different set of target moments may lead to a different set of calibrated parameters.

A poor fit does not diagnose its own cause. When a moment is badly matched, it is not obvious what the source of the poor fit is: a genuinely misspecified causal structure, weak identification (several very different parameterisations fitting about equally well), and sampling noise from a small or structurally unstable sample are all consistent with the same symptom. Telling them apart takes economic judgement beyond what the fitted moments alone can offer.

Summary and discussion

The Method of Simulated Moments (MSM) is a relatively easy-to-implement numerical method to calibrate the parameters of a model to empirical data. We have seen that calibrating even a small, well-known model raises difficult practical choices at every step: which empirical moments meaningfully identify which parameters; how to add an appropriate stochastic structure to the model; how to weight moments of very different precision and scale; and how to check that a numerical optimiser has actually found a good answer rather than a nearby one. Every one of these choices requires experimentation combined with careful reasoning. Unfortunately, many of the underlying questions do not have an obvious or objectively correct answer and require making judgement calls. This is why transparency and reproducibility are important: documenting and explaining each step carefully, so that other researchers can scrutinise the results and build on them. Below we give a simple recipe for implementing the method in practice – bearing in mind that some flexibility and independent reasoning is required at each step.

A key conclusion from our illustrative calibration of the simple 3-equation model against moments from US data over the period 1996 to 2007 is that the calibrated model does a mixed job at producing meaningful results. Put bluntly, it predicts a quasi permanent response to a one-off inflationary cost-push shock that the central bank is neither able nor willing to fight. While the 2021-23 inflation episode did involve a sustained increase in inflation, it also prompted a relatively aggressive response by many central banks. To what extent it was this monetary policy response rather than a change in the underlying drivers of inflation that brought inflation down is open to debate – however, inflation largely did stabilise again.

Unfortunately, the calibration exercise as such is silent on the question why the calibrated model generates a relatively implausible result. The culprit could lie in the way the model was calibrated or the model itself, or both. For example, a different data set, a different set of target moments, a different weighting matrix, and possibly even a different minimisation algorithm might yield a very different, possibly better, parametrisation of the model. At the same time, a different specification for the stochastic elements of the model, a different lag structure in the model (e.g. in the IS and/or Phillips curve) or a different assumption about inflation expectation formation might yield a better fit with the data. While this is perhaps not an entirely satisfactory conclusion, it shows that calibrating a model yields a lot of food for thought for further refinements of the model.

A simple recipe for calibrating economic models with the MSM

  1. Add a minimal stochastic structure to the model.
  2. Choose the free parameters of the model to be calibrated (i.e. those not already pinned down by a normalisation or fixed on other grounds).
  3. Choose target moments that are informative about the free parameters (at least as many moments as free parameters, i.e. \(n_m \geq n_p\)).
  4. Compute the chosen moments in the empirical data as well as in simulated model data, averaging the latter over many Monte Carlo runs.
  5. Choose a weighting matrix for the loss function.
  6. Minimise the loss function numerically, e.g. using R’s optim() function, restricting the search to a region of the parameter space that is economically plausible.
  7. Assess the fit of the calibrated model along several dimensions (e.g. whether the model matches the sign of each moment, percent deviations, moment-missing-rates).
  8. Check the robustness of the numerical solution (e.g. by varying the starting values for the minimisation algorithm, sample selection, data transformations, target moments, …).
  9. Interrogate a poor fit or implausible parameter estimate: do they reflect problems with the data, methodology, or theoretical model?.

References

Carlin, Wendy, and David Soskice. 2014. Macroeconomics. Instititions, Instability, and the Financial System. Oxford University Press.
Duffie, Darrell, and Kenneth J. Singleton. 1993. “Simulated Moments Estimation of Markov Models of Asset Prices.” Econometrica 61 (4): 929–52. https://doi.org/10.2307/2951768.
Franke, Reiner. 2018. “Competitive Moment Matching of a New-Keynesian and an Old-Keynesian Model.” Journal of Economic Interaction and Coordination 13: 201–39. https://doi.org/10.1007/s11403-016-0181-0.
Franke, Reiner. 2022. “An Empirical Test of a Fundamental Harrod-Kaldor Business Cycle Model.” Structural Change and Economic Dynamics 60: 1–14. https://doi.org/10.1016/j.strueco.2021.11.001.
Franke, Reiner, and Frank Westerhoff. 2012. “Structural Stochastic Volatility in Asset Pricing Dynamics: Estimation and Model Contest.” Journal of Economic Dynamics & Control 36 (8): 1193–211. https://doi.org/10.1016/j.jedc.2011.10.004.
Hamilton, James D. 2018. “Why You Should Never Use the Hodrick-Prescott Filter.” The Review of Economics and Statistics 100 (5): 831–43. https://doi.org/10.1162/rest_a_00706.
McFadden, Daniel. 1989. “A Method of Simulated Moments for Estimation of Discrete Response Models Without Numerical Integration.” Econometrica 57 (5): 995–1026. https://doi.org/10.2307/1913621.
Reissl, Severin. 2020. “Minsky from the Bottom up – Formalising the Two-Price Model of Investment in a Simple Agent-Based Framework.” Journal of Economic Behavior & Organization 177: 109–42. https://doi.org/10.1016/j.jebo.2020.06.012.
Wildauer, Rafael, Karsten Kohler, Adam Aboobaker, and Alexander Guschanski. 2023. “Energy price shocks, conflict inflation, and income distribution in a three-sector model.” Energy Economics 127B. https://doi.org/10.1016/j.eneco.2023.106982.

  1. See Franke and Westerhoff (2012), Franke (2018), Franke (2022), Reissl (2020), and Wildauer et al. (2023) for some applications.↩︎

  2. Franke (2018) is particularly close to the current chapter, but considerably more ambitious: he uses the MSM to compare a richer version of the New Keynesian model to a sentiment-driven Old-Keynesian model of business optimism and pessimism, matching up to 78 moments.↩︎

  3. If \(W\) is a diagonal matrix with entries \(w_i\), Equation 24.1 simplifies to \(J(\theta) = \sum_{i=1}^{n_m}w_i\big[m_i(\theta)-m_i^{emp}\big]^2\), where \(i=1,\dots,n_m\) indexes the moments.↩︎

  4. Alternatively, we could use the CBO’s GDPPOT potential-GDP series, which is a model-based estimate.↩︎

  5. The snapshot is available here.↩︎

  6. We leave the exploration of \(W=\hat{\Omega}^{-1}\) as an exercise to the reader.↩︎

  7. We compute this as \(\text{MMR}_i = \frac{100 \cdot |m_i - m_i^{emp}|}{2 \cdot SE_i}\), where \(SE_i\) is the bootstrapped standard error of moment \(i\).↩︎

  8. Since the model is linear, the response to any other shock size is simply a rescaled version of this one - e.g. the response to a shock half this size is half as large.↩︎

  9. The very long adjustment period can also be shown more formally. In Section 19.6, we show that the model’s single eigenvalue is given by \(\lambda = 1 - a_1 a_2 a_3\). With the calibrated parameters, this evaluates to \(\lambda \approx 0.9999\), i.e. a near unit-root.↩︎