Understanding and Implementing the Simulated Method of Moments (SMM)

The Simulated Method of Moments (SMM) is a powerful estimation technique. It allows us to estimate parameters of a model by matching moments from simulated data to those from observed data.

SMM is particularly valuable in scenarios where:

1. The model lacks closed-form analytical solutions
2. The data has limitations, such as only being available in aggregate format
3. The likelihood function is intractable
4. The model involves latent variables or complex dynamics that make traditional estimation methods impractical

The basic steps of SMM:

• Define the model and its parameters
• Choose moments to match
• Calculate moments from simulated and observed data
• Define an objective function that measures the distance between simulated and observed moments, weighted by a certain matrix
• Minimize the objective function to estimate parameters.

Let’s see a super super super basic example in R. We’ll build from here.

library(stats)
​
# Step 1: Define the true parameters (unknown in real applications)
true_mean <- 5
true_sd <- 2
​
# Generate observed data
set.seed(123)
observed_data <- rnorm(1000, mean = true_mean, sd = true_sd)
​
# Step 2: Choose moments to match (mean and variance)
observed_moments <- c(mean(observed_data), var(observed_data))
​
# Step 3 & 4: Function to simulate data and calculate moments
simulate_moments <- function(params, n = 1000) {
  simulated_data <- rnorm(n, mean = params[1], sd = params[2])
  c(mean(simulated_data), var(simulated_data))
  }
​
# Step 5: Objective function
objective <- function(params) {
  sim_moments <- simulate_moments(params)
  sum((sim_moments - observed_moments)^2)
  }
​
# Step 6: Minimize the objective function
result <- optim(par = c(0, 1), fn = objective, method = "L-BFGS-B", lower = c(-Inf, 0.01), upper = c(Inf, Inf))
​
# Print results
cat("Estimated mean:", result$par[1], "\n")
cat("Estimated sd:", result$par[2], "\n")
cat("True mean:", true_mean, "\n")
cat("True sd:", true_sd, "\n")

In actual applications of the SMM, researchers use a two-step optimization with optimal weighting matrix. This is a summary of how the procedure looks like:

• Perform an initial estimation using an identity matrix as the weighting matrix.
• Use the initial estimates to compute an optimal weighting matrix.
• Re-estimate the parameters using this optimal weighting matrix.

Let’s draft an example code in R to make this more concrete.

library(MASS)  # for mvrnorm function
library(numDeriv)  # for gradient and hessian calculations

# True parameters (unknown in real applications)
true_params <- c(mean = 5, sd = 2)

# Generate observed data
set.seed(123)
observed_data <- rnorm(1000, mean = true_params[1], sd = true_params[2])

# Calculate observed moments (mean and variance)
observed_moments <- c(mean(observed_data), var(observed_data))

# Function to simulate moments
simulate_moments <- function(params, n = 1000) {
  simulated_data <- rnorm(n, mean = params[1], sd = params[2])
  c(mean(simulated_data), var(simulated_data))
}

# Objective function
objective <- function(params, W) {
  sim_moments <- simulate_moments(params)
  diff <- sim_moments - observed_moments
  t(diff) %*% W %*% diff
}

# Two-step SMM estimation
smm_two_step <- function(initial_params) {
  # Step 1: Initial estimation with identity matrix
  step1 <- optim(par = initial_params, fn = objective, W = diag(2),
                 method = "L-BFGS-B", lower = c(-Inf, 0.01))
  
  # Compute optimal weighting matrix
  S <- cov(replicate(1000, simulate_moments(step1$par) - observed_moments))
  W_optimal <- solve(S)
  
  # Step 2: Re-estimate with optimal weighting matrix
  step2 <- optim(par = step1$par, fn = objective, W = W_optimal,
                 method = "L-BFGS-B", lower = c(-Inf, 0.01))
  
  return(list(params = step2$par, W = W_optimal))
}

# Run two-step SMM
result <- smm_two_step(c(0, 1))
print(result$params)

To estimate standard errors, bootstraps are often a good idea. This involves resampling our observed data, re-estimating the parameters, and calculating the standard deviation of these estimates. Let’s see how that may look like:

# Bootstrapping function
bootstrap_smm <- function(n_bootstrap = 1000) {
  bootstrap_estimates <- matrix(NA, nrow = n_bootstrap, ncol = 2)
  
  for (i in 1:n_bootstrap) {
    # Resample observed data
    boot_data <- sample(observed_data, replace = TRUE)
    boot_moments <- c(mean(boot_data), var(boot_data))
    
    # Define bootstrap objective function
    boot_objective <- function(params, W) {
      sim_moments <- simulate_moments(params)
      diff <- sim_moments - boot_moments
      t(diff) %*% W %*% diff
    }
    
    # Run SMM on bootstrap sample
    boot_result <- optim(par = result$params, fn = boot_objective, W = result$W,
                         method = "L-BFGS-B", lower = c(-Inf, 0.01))
    
    bootstrap_estimates[i,] <- boot_result$par
  }
  
  # Calculate standard errors
  se <- apply(bootstrap_estimates, 2, sd)
  return(se)
}

# Run bootstrap
bootstrap_se <- bootstrap_smm()

# Print results
cat("Estimated parameters:\n")
print(result$params)
cat("\nBootstrap standard errors:\n")
print(bootstrap_se)
cat("\nTrue parameters:\n")
print(true_params)

Now, some practical considerations.

1. Computational Intensity: Both the two-step procedure and bootstrapping can be computationally intensive. For large models or datasets, you might need to use parallel processing to speed up calculations - I recommend to always use parallel processing for bootstrap procedures. In many real-world applications, a standard laptop or desktop computer may not suffice. Researchers often need to utilize high-performance computing resources such as computing clusters or cloud computing services.
2. Initial Values: The choice of initial values can affect convergence. In practice, you might want to try multiple starting points. Sometimes you may need to impose certain boundaries - for example, a parameter estimate that is supposed to correspond to a standard deviation should be positive.
3. Moment Selection: The choice of moments can significantly impact estimation efficiency. In this example, I used mean and variance, but in more complex models, choosing informative moments is crucial. For identification, one needs at least as many moments as the number of parameters to estimate. In practice, it is a good rule of thumb to have a 2-3 to 1 ratio.
4. Weighting Matrix: Sometimes you may need to use a bootstrap procedure to estimate the inverse covariance-variance matrix as the optimal weighting matrix.

Remember, the key to successful SMM implementation lies in careful model specification, thoughtful moment selection, and robust optimization techniques. Happy simulation!