randy
randy

Reputation: 1081

R: quickly simulate unbalanced panel with variable that depends on lagged values of itself

I am trying to simulate monthly panels of data where one variable depends on lagged values of that variable in R. My solution is extremely slow. I need around 1000 samples of 2545 individuals, each of whom is observed monthly over many years, but the first sample took my computer 8.5 hours to construct. How can I make this faster?

I start by creating an unbalanced panel of people with different birth dates, monthly ages, and variables xbsmall and error that will be compared to determine the Outcome. All of the code in the first block is just data setup.

# Setup:
library(plyr)

# Would like to have 2545 people (nPerson). 
#Instead use 4 for testing.
nPerson = 4
# Minimum and maximum possible ages and birth dates
AgeMin = 10
AgeMax = 50
BornMin = 1950
BornMax = 1963

# Person-specific characteristics
ind = 
  data.frame(
    id = 1:nPerson,
    BornYear = floor(runif(length(1:nPerson), min=BornMin, max=BornMax+1)),
    BornMonth = ceiling(runif(length(1:nPerson), min=0, max=12))
  )

# Make an unbalanced panel of people over age 10 up to year 1986
# panel = ddply(ind, ~id, transform, AgeMonths = BornMonth)
panel = ddply(ind, ~id, transform, AgeMonths = (AgeMin*12):((1986-BornYear)*12 + 12-BornMonth))

# Set up some random variables to approximate the data generating process
panel$xbsmall = rnorm(dim(panel)[1], mean=-.3, sd=.45)
# Standard normal error for probit
panel$error = rnorm(dim(panel)[1])

# Placeholders
panel$xb = rep(0, dim(panel)[1])
panel$Outcome = rep(0, dim(panel)[1])

Now that we have data, here is the part that is slow (around a second on my computer for only 4 observations but hours for thousands of observations). Each month, a person gets two draws (xbsmall and error) from two different normal distributions (these were done above), and Outcome == 1 if xbsmall > error. However, if Outcome equals 1 in the previous month, then Outcome in the current month equals 1 if xbsmall + 4.47 > error. I use xb = xbsmall+4.47 in the code below (xb is the "linear predictor" in a probit model). I ignore the first month for each person for simplicity. For your information, this is simulating a probit DGP (but that is not necessary to know to solve the problem of computation speed).

# Outcome == 1 if and only if xb > -error
# The hard part: xb includes information about the previous month's outcome
start_time = Sys.time()
for(i in 1:nPerson){
  # Determine the range of monthly ages to loop over for this person
  AgeMonthMin = min(panel$AgeMonths[panel$id==i], na.rm=T)
  AgeMonthMax = max(panel$AgeMonths[panel$id==i], na.rm=T)
  # Loop over the monthly ages for this person and determine the outcome
  for(t in (AgeMonthMin+1):AgeMonthMax){
    # Indicator for whether Outcome was 1 last month
    panel$Outcome1LastMonth[panel$id==i & panel$AgeMonths==t] = panel$Outcome[panel$id==i & panel$AgeMonths==t-1] 
    # xb = xbsmall + 4.47 if Outcome was 1 last month
    # Otherwise, xb = xbsmall
    panel$xb[panel$id==i & panel$AgeMonths==t] = with(panel[panel$id==i & panel$AgeMonths==t,], xbsmall + 4.47*Outcome1LastMonth)
    # Outcome == 1 if xb > 0
    panel$Outcome[panel$id==i & panel$AgeMonths==t] =
      ifelse(panel$xb[panel$id==i & panel$AgeMonths==t] > - panel$error[panel$id==i & panel$AgeMonths==t], 1, 0)
    }
  }
end_time = Sys.time()
end_time - start_time

My thoughts for reducing computer time:

Upvotes: 1

Views: 813

Answers (1)

Weihuang Wong
Weihuang Wong

Reputation: 13108

One approach is to use a split-apply-combine method. I take out the for(t in (AgeMonthMin+1):AgeMonthMax) loop and put the contents in a function:

generate_outcome <- function(x) {
  AgeMonthMin <- min(x$AgeMonths, na.rm = TRUE)
  AgeMonthMax <- max(x$AgeMonths, na.rm = TRUE)
  for (i in 2:(AgeMonthMax - AgeMonthMin + 1)){
    x$xb[i] <- x$xbsmall[i] + 4.47 * x$Outcome[i - 1]
    x$Outcome[i] <- ifelse(x$xb[i] > - x$error[i], 1, 0)
  }
  x
} 

where x is a dataframe for one person. This allows us to simplify the panel$id==i & panel$AgeMonths==t construct. Now we can just do

out <- lapply(split(panel, panel$id), generate_outcome)
out <- do.call(rbind, out)

and all.equal(panel$Outcome, out$Outcome) returns TRUE. Computing 100 persons took 1.8 seconds using this method, compared to 1.5 minutes in the original code.

Upvotes: 1

Related Questions